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ABSTRACT 


The  linear  theory  of  viscoelasticity  is  used  to 
study  the  effects  of  anelasticity  on  seismic  body  waves 
propagating  through  a  layered  homogeneous  isotropic  medium. 
The  mathematical  formulation  describing  the  propagation  and 
attenuation  of  plane  waves  in  a  homogeneous  isotropic  linear 
viscoelastic  medium  is  given.  The  particle  motion  associa¬ 
ted  with  plane  viscoelastic  waves  is  examined,  and  is  seen 
to  be  elliptical  for  P-SV  waves  and  linear  for  SH  waves. 

A  conservation  of  energy  equation  for  steady-state  harmonic 
viscoelastic  waves  is  derived  from  first  principles  and  is 
used  to  calculate  various  time-averaged  energy-related 
quantities  such  as  energy  fluxes,  potential  and  kinetic 
energy  densities,  and  dissipation  rates.  Reflection  and 
transmission  coefficients  for  plane  waves  impinging  upon  a 
plane  boundary  separating  two  linear  viscoelastic  media 
are  calculated.  The  coefficients  for  SH  waves  are  plotted 
and  compared  with  the  corresponding  coefficients  for  the 
perfectly  elastic  case.  Asymptotic  ray  theory  is  extended 
to  homogeneous  isotropic  linear  viscoelastic  media,  and  the 
geometrical  spreading  factor  for  a  given  viscoelastic  ray 
generated  by  a  surface  point  source  is  calculated.  Synthe¬ 
tic  seismograms  for  SH  body  waves  are  computed  for  a  plane¬ 
layered  crustal  model  in  both  the  elastic  and  anelastic 
cases,  using  a  ray  theory  approach.  Both  the  teleseismic 
case,  and  the  case  of  a  point  source  at  the  surface  are 
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treated . 


The  seismograms  for  the  anelastic  medium  exhibit 
attenuation  and  waveform  spreading. 


amplitude 
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CHAPTER  1 


INTRODUCTION 


It  is  becoming  increasingly  important  to  include  the 
effects  of  anelasticity  in  the  computation  of  synthetic 
seismograms  and  in  other  seismic  modelling  techniques. 
Anelastic  effects  such  as  the  attenuation  of  seismic  waves , 
are  often  modelled  by  the  use  of  empirical  laws  deduced 
from  seismic  data.  However,  to  facilitate  the  accurate 
interpretation  of  seismic  field  records,  and  to  better 
understand  the  behaviour  of  seismic  waves  propagating  in  an 
anelastic  medium,  it  is  desirable  to  treat  the  anelasticity 
from  a  more  physical  viewpoint. 

In  this  regard,  a  possible  approach  may  be  to  use 
the  linear  theory  of  viscoelasticity  to  theoretically 
model  the  anelasticity.  In  this  theory,  an  anelastic 
material  is  said  to  have  "memory".  In  other  words,  the 
state  of  the  material  at  a  given  time  depends  on  the  history 
of  the  material  up  to  the  given  time.  For  instance,  in  the 
one-dimensional  case,  the  stress  on  a  solid  at  a  given  time 
t  depends  on  all  the  strain  forces  that  have  acted  on  the 
solid  in  the  past,  up  to  the  time  t.  In  mathematical  terms. 


a  (t) 


—  00 


r (t-T) 


de  ( x ) 
dr 


dr 


where  a  is  the  stress,  e  is  the  strain,  and  r  is  known  as 
the  relaxation  function,  and  is  a  characteristic  function 
of  the  solid.  This  equation  is  known  as  the  Boltzman 
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after-effect  equation.  Various  spring-dashpot  models  can 
be  constructed  to  describe  different  types  of  viscoelastic 
solids.  For  each  different  model  the  relaxation  function 
can  be  calculated,  and  the  loss  factors,  which  are  measures 
of  internal  friction,  can  be  calculated  from  the  relaxation 
function.  Detailed  expositions  of  the  linear  theory  of 
viscoelasticity  can  be  found  in  such  books  as  Fung  (1965), 
Christensen  (1971)  ,  Lockett  (1972)  and  Fliigge  (1975)  . 

The  main  purpose  of  this  thesis  is  to  examine  the 
nature  of  seismic  waves  propagating  in  anelastic  media  and 
to  investigate  and  develop  methods  of  computing  synthetic 
seismograms  for  plane  layered  anelastic  media  using  a  ray 
theory  approach.  The  linear  theory  of  viscoelasticity  is 
used  to  model  the  anelasticity .  The  media  are  assumed  to 
be  homogeneous  and  isotropic.  The  numerical  values  of  the 
loss  factors  (i.e.,  the  reciprocal  Q's)  for  the  various 
layers  are  assumed  to  be  known,  and  are  taken  as  input  para 
meters  in  the  calculations.  The  main  advantage  in  this  is 
that  all  calculations  are  then  independent  of  any  viscoelas 
tic  spring-dashpot  model. 

Since  a  ray  theory  approach  is  used,  it  is  necessary 
to  compute  reflection  and  transmission  coefficients  for 
plane  waves  impinging  upon  a  plane  boundary  separating  two 
anelastic  media.  Some  interesting  aspects  of  the  behaviour 
of  plane  waves  propagating  in  layered  anelastic  materials 
are  revealed  by  the  examination  of  these  coefficients  and 
of  the  reflected  and  transmitted  waves. 
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Synthetic  seismograms  are  computed  for  both  the 
teleseismic  case,  and  the  case  of  a  point  source  at  the 
surface.  In  the  latter  case,  asymptotic  ray  theory  is 
extended  to  include  linear  viscoelasticity  and  used  to  cal¬ 
culate  the  geometrical  spreading  factor  for  a  given  ray  in 
a  layered  anelastic  medium. 

The  results  of  the  viscoelastic  calculations  are 
often  compared  with  the  analogous  results  for  the  perfectly 
elastic  case,  in  order  to  better  understand  the  effects  of 
anelasticitv .  The  elastic  results  are  easily  obtained 
since  the  theory  of  elasticity  is  contained  within  the 
theory  of  viscoelasticity  as  the  special  case  of  no  dissi¬ 
pation  . 


CHAPTER  2 


PLANE  WAVE  PROPAGATION  AND  ATTENUATION 


In  this  chapter,  we  give  the  mathematical  formulation 
necessary  to  describe  the  propagation  and  attenuation  of 
plane  seismic  waves  in  a  homogeneous  isotropic  linear  visco¬ 
elastic  medium  (see  Borcherdt  (1973,  1977),  Silva  (1976) 
and  Buchen  (1971a) ) . 

The  stress-strain  relation  for  a  homogeneous  iso¬ 
tropic  linear  viscoelastic  medium  is  given  by 


a  .  .  ( t) 
1 D 


t 

c 


ID 


X(t-T) 


dekk(T) 

dx 


dx  +  2 


y (t-x) 


de^  .  (x) 
dx 


dx 


=  6.  .A(t)  *  de,,  (t)  +  2y(t) 
1  3  kk 


*  de . . ( t ) 
iD 


(2.1) 


where  the  symbol  *  denotes  the  Stieltjes  convolution  (see, 
for  example,  Christensen  (1971)  and  Fung  (1965)).  X  and 
y  are  the  Lame  parameters,  and  e^  is  the  strain  tensor 
given  by 


e .  .  =  T  ( u .  .  +  u .  .  ) 
ID  2  1,3  D  / 1 


(2.2) 


where  u^,  i  =  1,2,3  is  the  displacement  vector.  The  Ein¬ 
stein  summation  convention  is  used  in  (2.1)  ,  and  in  the 
equations  that  follow.  As  in  the  case  of  perfect  elasti¬ 
city,  the  viscoelastic  stress-strain  relation  can  be  split 
into  bulk  and  shear  components  as  follows : 
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tfkk(t)  =  3K(t)  *  dekk(t) 

c^i j  ( t )  =  2y(t)  *  de^j  (t)  ,  i  f  j 

where 

K(t)  =  A  ( t)  +  |y(t) 

Substituting  (2.1)  into  the  equation  of  motion 


(2.3) 


(2.4) 


a 


ij  ,  j 


P 


u . 
i 


(2.5) 


yields,  after  some  calculation. 


[A (t) +y (t) ]*d(V (V-u))  +  y(t)*d(V2  u)  =  p  u  (2.6) 


in  Cartesian  coordinates.  Taking  the  Fourier  transform 
(see  Appendix  1) ,  we  get 


(A+M)V(V*u)  +  M  V  u  =  -  poj  u 


(2.7) 


where 

-> 

u 

A 

M 


-*  -ioot 
u  e  dt 


A(co)  =  iw 


M(oo)  =  im 


o 

CO 


J0 


A ( t )  e  dt 


v  -iwt 
y  ( t)  e  dt 


(2.8) 


n 
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,  ->■  t 

We  can  write  u  in  terms  of  Helmholtz  potentials  as 


u  =  V(f)  +  V  x  ip 


(2.9) 


and  inserting  this  in  (2.7)  gives 


2  - 

V  <J> 


.  2  - 
+  kp  (J) 


=  0 


V 


2 


ip  =  0 


where 


(2.10) 


kp2  =  co2/a2  ,  a2  =  (A+2M)/p 

2  2  2  2 
ksz  =  w /e  ,  r  =  m/p 


(2.11) 


and  where  p  is  the 
velocities  a  and  3 
The  general 


density.  We  can  see  that  the  P  and  S 
are  now  complex,  and  frequency-dependent. 
Helmholtz  equation  can  be  written  as 


2+  2  + 

V  F  +  k  F  =  0 


(2.12) 


which  applies  to  both  P  waves  (k  =  k  ,  F  =  <j>)  and  S  waves 
(k  =  k  ,  F  =  \p)  .  Its  solution  is  a  general  plane  wave 

O 

given  by 


F 


-A«r  -iP*r  ^  -ik*r 
F  e  e  =  F  e 

o  o 


(2.13) 


* 
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where 


k  =  P  -  iA.  (2.14) 

P  is  the  propagation  vector,  and  is  perpendicular  to  the 
planes  of  constant  phase  defined  by  P*r  =  constant.  A  is 
the  attenuation  vector  and  is  perpendicular  to  the  planes 
of  constant  amplitude  defined  by  A*?  =  constant.  The 
velocity  of  the  planes  of  constant  phase  is  given  by 

LO  * 

v  =  P  (2.15) 

I  p  I 

^  .  ,  ,  , 
where  P  is  a  unit  vector  in  the  directon  of  P. 

From  (2.14)  we  have 

2  -*  2  2 

k  =  k*k  =  P  -  A  -  2i  P-A  (2.16) 


where 


P  «A 


PA  cos  y 


(2.17) 


where  y  is  the  angle  between  P  and  A,  and  P  and  A  are  the 
t 

magnitudes  of  P  and  A.  If  P  and  A  are  parallel  (y  =  0) , 
the  wave  is  called  homogeneous.  If  P  and  A  are  not  parallel 
(y^O) ,  the  wave  is  called  inhomogeneous.  To  ensure  that  the 
amplitude  of  the  wave  does  not  increase  in  the  direction  of 
propagation,  we  must  have,  in  general. 


. 
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o  *  M  *  £ 


(2.18) 


Solving  (2.16)  and  (2.17)  for  P  and  A,  we  get 


( 


P2  =  i-le(k2)  +  { [Re  (k2)  ] 2  +  R-V-MJ  .  }* 
\  cos  y 


AZ  =  j[  -Re  ( k  2 )  +  {  [  Re  ( k  2 )  ]  2  +  -^-im  (k  J-L . 

cos  Y  j 


(2.19) 


We  can  also  express  P  and  A  in  a  more  amenable  form  in  terms 
of  the  loss  factors  and  homogeneous  phase  speeds  of  the 
medium,  which  are  assumed  to  be  known,  as  shown  below. 

Following  Borcherdt  (1973),  we  define  the  loss 
factor  for  P  and  S  waves  as 

Qs  =  Im(M)/Re(M)r  S  waves 

Q_1  =  (2.20) 

Qp  1  =  Im (A+2M) /Re ( A+2M) ,  P  waves. 

2  .  -1 
k  can  then  be  written  m  terms  of  Q  as 


2 ( 1-iQ  1)  ] 

l  +  /l+Q~2  J 


(2.21) 


for  both  P  and  S  waves  (i.e.,  k=kg,  Q=Qs'v^=vt^s  ^or  s 
waves,  and  k=kp ,  Q=Qp,  VH=VHP  for  p  waves),  where  vR  is  the 
homogeneous  phase  speed  of  the  wave.  Substituting  (2.21) 
into  (2.19)  yields 


■ 


9 


2 

CO 

/ 

1  + 

/  -2  2 
/l+Q  /cos  y 

2 

V 

H 

1 

+  /l+Q-2 

2 

CO 

-1 

+  /l+Q  ^/cos^y 

v  2 

H 

\  1 

+  /l+Q-2 

(2.22) 


for  both  P  and  S  waves . 

2 

If  the  medium  is  non-dissipative  (Im(k  )  =  0)  , 
equation  (2.17)  shows  that  either  A  =  0  or  y  =  ±  tt/2. 
Conversely,  if  A  -  0  or  y  =  ±  tt/2  the  medium  is  non-dissi¬ 
pative.  Hence,  we  can  say  that  the  medium  is  dissipative 
if  and  only  if  A  /  0  and  y  /  ±  tt/2  (i.e.,  |y|  <  tt/2). 

This  shows  that  the  only  type  of  inhomogeneous  plane  wave 
that  propagates  in  a  non-dissipative  medium  (i.e.,  the  wave 
with  A  /  0  and  y  =  ±  tt/2)  does  not  propagate  in  a  dissipa¬ 
tive  medium,  and  vice  versa.  This  is  a  basic  difference 
between  the  theories  of  elasticity  and  viscoelasticity. 

Using  the  subscripts  "IH"  for  "inhomogeneous"  and 
"H"  for  "homogeneous",  we  can  see  from  equations  (2.15)  and 
(2.19)  or  (2.22)  that 


i  i 

i  -*■  i 

V„„ 

<  vl 

1  IH  1 

1  H  1 

i  p 

i  ->  i 

A 

>  A 

1  "IH  1 

1  H 1 

(2.23) 


i.e.,  the  phase  speed  of  an  inhomogeneous  wave  is  less  than 
that  of  a  homogeneous  wave,  and  the  maximum  attenuation  of 
an  inhomogeneous  wave  is  greater  than  that  of  a  homogeneous 


■ 


- 


wave . 


P  and  A  can  be  calculated  uniquely  from  the  medium 
parameters ,  for  a  specified  y,  from  (2.19)  or  (2.22). 

However,  the  velocity  and  attenuation  of  an  inhomogeneous 
wave  in  a  non-dissipative  medium  (A  /  0,  y  =  ±  tt/2)  cannot 
be  uniquely  determined  from  (2.19)  or  (2.22).  They  can 
only  be  determined  once  the  boundary  conditions  have  been 
specified.  For  an  example  of  this  non-uniqueness  in  the 
case  of  elasticity,  consider  the  following  two  waves:  an 
interface  P-wave  generated  by  an  SV  wave  incident  on  a 
plane  interface  at  an  angle  greater  than  the  P-critical 
angle,  and  a  Rayleigh  wave  on  an  elastic  half-space.  Both 
waves  have  y  =  tt/2,  but  their  phase  velocities  are  different. 


CHAPTER  3 


PARTICLE  MOTION  AND  ENERGY 

3 . 1  Particle  Motion 

Borcherdt  (1973)  and  Buchen  (1971a)  have  shown  that 
the  particle  motion  for  P-SV  steady-state  harmonic  plane 
viscoelastic  waves  is  elliptical,  rather  than  linear  as  in 
the  case  of  perfect  elasticity. 

This  is  demonstrated  as  follows.  First,  let  us 
employ  the  usual  steady-state  harmonic  condition  on  the 
displacement , 


and  insert  this  into  the  time-domain  equation  of  motion 
(2.6).  We  obtain 


(A4M) V (V«u) 


+  MV2u 


(3.2) 


where  A  and  M  are  given  by  (2.8)  .  Writing  u  in  terms  of 
Helmholtz  potentials 


u  =  V(f>  +  Vx\ jj 


and  inserting  this  into  (3.2)  yields 


(3.3) 


11 


' 


12 


V24>  +  kp2ct>  =  0 

2-y  2~> 

V  ip  +  k  ip  =  0 

O 


(3.4) 


where  k  and  k  are  given  by  (2.11)  . 

±r  iD 

For  P  waves,  the  solution  for  (p  is  given  by 


B 


i  (cot-k  •  r) 
e 


(3.5) 


where  B  is  a  complex  constant,  k  =  P-iA,  and  where  the 
subscript  "P"  has  been  suppressed  in  all  the  parameters. 
The  displacement  field  is  then  given  by 


u  =  V(J> 


=-i&B 


i (wt-k • r ) 


=  -ie  A’r  (kB) 


k  i (mt-P • r ) 
k  e 


=  C 


k  ic (t) 
k 


(3.6) 


where 


C 


kB 


-> 

•  r 


^(t)  =  oot  -  ?«r  +  arg(kB)  -  tt/2 


The  quantity  k/k  can  be  written  as 


(3.7) 


'  h  ' 
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k 

k 


->  .  -> 


P  -  lA 
kR+ikI 


(3.8) 


where 


k?-ki 

|  =  _* _ L- 


ki?+kRS 


(3.9) 


The  subscripts  "R"  and  "I"  refer  to  real  and  imaginary 

parts.  Inserting  these  in  (3.6),  and  taking  the  real  part 

— f 
of  u  to  get  the  actual  physical  motion,  yields 


uR  =  C  [  f1  cos  c(t)  +  %2  sin  C(t)] 


(3.10) 


Also,  using  (2.16)  and  some  rules  of  complex  algebra,  it  is 
easy  to  show  that 
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1 


K 


(3.11) 


Hence  and  f2  are  perpendicular,  and  ?]_  >  ~  °* 

Setting 


x 


-> 

u 


1 


V^2 

C?2 


5,  cos  c (t) 


(3.12) 


Y 


£2  sin  C(t) 
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reveals  that 


which  is  an  ellipse  with  major  axis  £  ,  and  minor  axis 
The  motion  is  in  the  plane  of  P  and  A,  and  the  direction 
of  rotation  is  from  P  to  A.  For  a  homogeneous  P  wave 
(y  =  0)  ,  it  is  easy  to  show  that  =  1  and  =  0,  hence, 
the  elliptical  particle  motion  degenerates  to  straight-line 
or  linear  motion  of  a  longitudinal  nature,  i.e.  parallel 
to  the  direction  of  propagation  P. 

A  similar  calculation  can  be  performed  for  the  SV 
wave,  with  similar  results.  The  solution  to  (3.4)  for  ip 
representing  the  SV  case  is  given  by 


-y  -> . 


-y 

B 


i (wt-k • r ) 


(3.14) 


where  B  =  Bn  with  B  a  complex  constant  and  n  a  real  unit 
vector,  and  where  the  subscript  "S"  has  been  suppressed  in 
all  the  parameters.  The  subsidiary  condition 

V  •  $  =  0  (3.15) 

implies  that 


->  -y  ,  -y  ,  -y .  -y  _ 
k  •  B  =  (P  -  lA)  •  B  =  0 


(3.16) 
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Hence, n  is  perpendicular  to  P  and  A.  The  displacement 
field  is  given  by 


u  =  V  x 


=  -i  (k  X  S)  ei(“t-k-r) 


.  -A  •  r 

-ie  ‘  (kB) 


it 

kxn 


i (mt-P *r) 


r  k 

C  7-  x  n 

k 


iC  (t) 


(3.17) 


where 


C 


kB 


e 


-A* 


-> 

r 


C(t)  =  rnt  -  P*r  +  arg(kB)  -  tt/2 


(3.18) 


Following  the  same  procedure  as  for  P  waves,  using 


u  =  Re  (V  x  $) 

K, 


(3.19) 


we  again  obtain  equations  (3.10)  -  (3.13),  where  for  SV 
waves , 


k 

—  X 

k 


n 


-y 

=  5- 


(3.20) 


with 
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-> 


(kRP-k].A)x  n 


fk  P+k  A)  x  n 


(3.21) 


For  homogeneous  waves,  £  =  1  and  =  0 ,  so  the  particle 

motion  degenerates  to  linear  transverse  motion,  perpendi- 

— y 

cular  to  P. 

Anelastic  P-SV  waves  exhibit  elliptical  particle 
motion.  However,  for  anelastic  SH  waves,  Borcherdt  (1977) 
has  shown  that  the  particle  motion  is  linear  in  the  direc- 
tion  perpendicular  to  the  plane  containing  P  and  A,  i.e., 
the  x 2  or  y  direction,  if  we  assume  P  and  A  lie  m  the 
Xl-X3  or  X-Z  P^ane *  This  is  the  same  motion  as  in  the  case 
of  perfect  elasticity.  The  solution  to  (3.4)  for  i  repre¬ 
senting  the  SH  case  is  given  by 


->  ->■ 
ijj  =  B 


i (wt-k*r) 
e 


(3.22) 


where 

B  =  x^  +  x^  (3.23) 

where  z^,z^  are  arbitrary  complex  numbers,  and  x^,x^  are 
real  unit  vectors  in  the  x^  and  x^  directions.  Calculating 
the  displacement  field,  we  obtain 


1 1 
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->  -> 
u  =  V  x  \p 


=  -  l 


(kxB) 


i  (wt-k • r) 


_  i  (ajt-k*  r) 
-  D  e 


/\ 


(3.24) 


where 


D  =  ( z ~  zi^3)  *  &  +  iP)  (3.25) 

Taking  the  real  part  of  u  to  obtain  the  actual  physical 
motion,  yields 

uR  =  |d  !  e  r  cos  (cot-P*r  +  arg  (d  )  )  (3.26) 

which  represents  a  linear  motion  in  the  direction. 

3 . 2  Energy 

In  this  section,  we  calculate  various  energy-related 
quantities  such  as  kinetic  and  potential  energy  densities, 
dissipation  rates,  and  energy  fluxes  for  steady-state 
harmonic  plane  waves  propagating  in  a  linear  viscoelastic 
medium.  Such  energy  calculations  have  been  carried  out 
previously  by  Borcherdt  (1973)  and  Buchen  (1971a),  who  used 
two  different  approaches.  The  calculations  in  this  section 
combine  the  general  aspects  of  these  two  papers  to  produce 


1 


I  *  % 
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a  more  general  approach  to  the  calculation  of  energy- 
related  quantities.  The  approach  is  also  different  from 
the  above-mentioned  papers  in  that  the  equation  for  the 
conservation  of  energy  is  calculated  from  first  principles . 

The  stress  tensor  for  the  steady-state  harmonic 
case  is  obtained  by  substituting  (3.1)  into  (2.1) .  The 
result  is 

a.  .  =  A(ai)  0  6.  .  +  2  M  ( oo )  e.  .  (3.27) 

ID  13  13 


where 


0  =  e.  .  =  V  •  u 

kk 


(3.28) 


To  obtain  an  equation  for  the  conservation  of  energy, 
we  begin  by  considering  the  power  P,  or  rate  of  work  done 
on  a  volume  V  with  boundary  S.  Neglecting  body  forces,  P 
is  given  by 


P 


dS 


S 


(3.29) 


where  t  is  the  traction  and  v  is  the  velocity.  The  trac¬ 
tion  is  given  by 


t  = 


-> 

-> 


a 


-> 

•  v 


(3.30) 


where  a  is  the  stress  dyad  and  v  is  a  unit  vector  normal  to 


*3 
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the  boundary.  Using  the  symmetry  of  the  stress 


a  .  .  =  a  .  . 
ID  D 1 


and  the  dyad-vector  rules 


(X  •  Y) 


£  X.  .  Y  . 
D  1D  D 


(Y  •  X)  ± 


£  Y  .  X  .  . 
D  : 


and  Gauss's  theorem,  v/e  have 


V  = 


_  f 


• 

( a  •  v )  •uD  dS 
j  R  R 

S 


a .  v .  u.  dS 
ljR  j  lR 


(a^uR)  -v  dS 


-V  • 

V*  (  aR*uR)  dV 


V 


(0ijR  UiR>'jdV 


V 


^•aij,jRUiR  +  aijR  Ui,jR^  dV 


V 


Using  the  equation  of  motion  (2.5) ,  and  the  easi 
strated  fact 


°i j R  Ui,jR  °ijR  SijR 


tensor ,  i . e . , 

(3.31) 

(3.32) 


(3.33) 

ly  demon- 


(3.34) 
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equation  (3.33)  becomes 


P  = 


[p“iR  UiR  +  °ijR  hjR]  dV 


V 


_3_ 
3 1 


1  *  * 

^pu.^u.ra  dV  + 
2  lR  lR 


V 


a .  e . .  dV 
13R  13R 


V 


(3.35) 


Using  (3.33),  we  then  obtain 


3 1 


2PUiRUiR 


"1 

dV  1  + 


a . .  e . .  dV 
13R  ijR 


• 

-*  -> 

'  °R  UR 


V 


dS 


(3.36) 


where 


1  •  • 

75-pu.  u.  =  kinetic  energy  density  =  KED 
Z.  1 K  1 K 


a.  .  _e.  . _  =  rate  of  change  of  strain  energy 
13R  ljR  ^  ^-L 


density 


(3.37) 


->  -> 


-o_  •  u_  =  energy  flux  vector,  or  intensity  =  I 

K  K 


Using  (  3 . 2  7  > ,  (3.32)  and  the  dyad  rule 


( Vu)  .  .  =  u  .  =  3  .  u  .  =  u  .  . 

13  3xi  3  13  3  / 1 


(3.38) 


the  rate  of  change  of  strain  energy  density,  and  the 
intensity  I  can  be  expanded  and  expressed  in  terms  of  A,  M 
and  the  displacements.  When  this  is  done,  (3.36)  finally 


1 
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becomes,  after  some  algebra 


_9_ 
9 1 


(KED  +  PED)  dV 


+  V  dV 


i-v  dS  (3.39) 


V 


V 


This  is  the  equation  expressing  the  conservation  of  energy 
for  steady-state  harmonic  viscoelastic  waves.  KED  is  the 
kinetic  energy  density  and  PED  is  the  potential  energy 
density.  V  is  the  rate  of  energy  dissipation  per  unit 

~y  -}■  ->■  , 

volume.  I  =  1^  +  is  the  energy  flux,  where  1^  is  the 
rate  at  which  the  material  inside  (outside)  the  surface  S 


is  doing  work  on  the  material  outside  (inside),  and  is 


the  rate  at  which  energy  is  convected  across  the  surface  S. 
The  mathematical  expressions  for  these  terms  are 


v  -  l  [aA2  +  2MI 


(3.40) 


I 


-v 

I 


2 


h  [AI®R5R  +  MI  (SR  •  v3R  +  V^R  '  “r>] 


■ 


. 
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We  now  calculate  the  time  averages  of  several 
quantities  for  P,  SV  and  SH  waves,  making  use  of  the  well 
known  facts 

<  cos2  [cot  -  f(r)]  >  =  1/2 

<  sin2  [cot  -  f  (r)  ]  >  =  1/2  (3.41) 

<  sin  [cot  -  f  (r)  ]  cos  [cot  -  f  (r)  ]  >  =  0 


where  <  ...  >  denotes  the  time  average  over  one  cycle. 


3.2.1  P  Waves 


The  scalar  potential  c p  for  P  waves  is  given  by 
(3.5) ,  where  k  =  P  -  lA.  The  displacement  is  then  given 
by 


u 


Vcf) 


-ilcB 


i  (cot-k  •  r ) 
e 


(3.42) 


From  (3.27),  we  obtain 


a .  . 

13 


(cot-k  *r) 


(3.43) 


where 


£.  .  =  Ak26 .  .  +  2Mk.k  . 
13  13  13 


(3.44) 


1 


If 
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The  subscript  i  (i=l,2,3)  is  not  to  be  confused  with  the 


complex  number  i.  The  mean  energy  flux  is  then 

given  by 

<1 .  >  =  -  <o . u.  > 
l  ljR  jR 

(3.45) 

Using  (3.41)  -  (3.43),  this  becomes,  after  some 

algebra , 

<I.>  =  4  03 1  B  !  2  e"2A*r  Re  [  C  •  •  k.*] 

1  2  1  1  ID  j 

(3.46) 

where  the  superscript  *  denotes  the  complex  conjugate. 
Using 


2  2  2 

Ak  =  pto  -  2Mk 

(3.47) 

C . .  can  be  written  as 
ID 


£.  .  =  poo26  .  .  +  2  M  ( k  .  k  .  -  k2  6  .  .) 

ID  ID  1  D  ID 

(3.48) 

Hence 


C.  .k.*  =  po)2k .  *  +  2M  (k.k.  -  k26.  .)  k.* 

ID  D  1  ID  ID  D 

(3.49) 

It  can  be  shown  that 


(k.k. -k26 .  .  )  k.*  =  2 (A . +iP  .  )  (P.A.-P.A.) 

'  1  j  ^  j  ^  “j  j  -j 

=  [2 (A  x  p)  x  (A  +  iP) ]± 

(3.50) 

(3.50) 
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Using  (3.49)  and  (3.50)  in  (3.46),  we  finally  obtain  for 
the  mean  energy  flux. 


<I> 


B 


2 


-2A  •  r 
e 


2=> 
poj  P 


+  4 (PxA) x (M  P-MrA) ] 


(3.51) 


For  homogeneous  waves  (PxA  =  0),  this  reduces  to 


<I> 


H 


2  -2A • r 

e 


-> 

P 


(3.52) 


and  for  a  perfectly  elastic  material  (Mj  =  0,  P*A  =0),  it 
reduces  to 


<I> 


B 


2 


-2A  •  r 
e 


[pto2+  4M  A2] 

K 


p 


(3.53) 


where  the  vector  identity 

,  ->  v  ->v  ->  .  .  “>■ 
(PxA)  x  A  =  (P  •  A)  A  -  ( A  •  A)  P 


(3.54) 


was  used.  From  (3.51)  -  (3.53)  we  see  that  the  direction 
of  the  mean  energy  flux  is  not,  in  general,  the  same  as  the 
direction  of  propagation,  except  for  anelastic  homogeneous 
waves,  and  elastic  waves. 

Next,  we  calculate  the  mean  kinetic  energy  density. 
Using  (3.41)  and  (3.42),  it  is  easily  calculated,  and  is 
given  by 
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<KED>  = 


<  2  PUiR  UiR> 


1  2  .  .2  -2A  •  r  ,_2  ,  _  2 . 

j  poo  |  B  |  e  (P  +A  ) 


(3.55) 


The  calculation  of  the  mean  potential  energy  density 
is  more  complicated.  Using  the  expression  given  for  PED 
in  (3.40),  and  using  (3.41)  and  (3.42),  we  obtain,  after 
some  algebra, 


<PED>  =  <  |  Ar6r2  +  Mr  e.jR  e. jR> 


_|2  -2A*r  r  1.  | .  2  |  2  1,,  ,2  2, 2-, 
B  |  e  [  jA^k  |  +  ^ R  (P  +A  )  ] 


(3.56) 


Equation  (3.47)  implies 


Ar  =  poo2  Re  (1/k2)  -  2Mr 


2  (k% 

=  pw  - -  -  2M 

i,2>2 


R 


(3.57) 


Hence , 


0,9  99  ,  9  ,  9 

Ar  \k~\Z  =  (k  ) R  -  2MR | k  | 


(3.58) 


2  ,2  ^2,  .,2.2 

=  poo  (P  -A  )  -  2MR|k  | 


and  so  <PED>  becomes 


I 
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<PED>  = 


1 

4 


B 


-2A  ■ 


-y 

r 


[poo2  (P2-A2)  +2M 


R 


{(p2+a2)2- 


k2|2)] 


(3.59) 


Using  (2.16)  ,  one  can  show  that 


2  2  2 
(P  +A  ) 


,  2  |  2  „ 
k  =4 


->■  -y  |  2 
P  x  A 


(3.60) 


Hence,  we  finally  obtain 

<PED>  =  j  |  B  |  2  e"2A*r  [  poo2  (P2-A2 )  +  8MR  |PxA|2]  (3.61) 

Note  that  <PED>  ^  <KED>  unless  A  =  0  which  occurs  only  for 
homogeneous  waves  in  elastic  media.  Using  (3.55)  and 
(3.61),  we  can  also  obtain  the  total  mean  energy  density 
<  E  >  : 


<  E  >  =  <KED>  +  <PED> 


1  ,,2  -2A • r 

2  |B|  6 


o  o 

[poo  P  + 


4MC | Pxa|  ] 

K 


(3.62) 


We  can  now  derive  an  interesting  relationship  between 
<E>  and  <I>.  Using  the  vector  identity 


-y  -y  ->  -y 

A* (BXC)  =  B 


(CxA)  = 


-y  -y  -y 

C • (AXB) 


(3.63) 


we  can  show  that 


H 


- 
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p.{(pxA)  x  (H  J-PJ)  }  =  M_  I  PxA  I  2  (3.64) 

-1  K  R 

Using  this  in  (3.51)  gives 

<E>  =  -  P  •  <I>  (3.65) 

to 

i.e.,  the  mean  total  energy  density  is  proportional  to  the 
component  of  the  mean  energy  flux  along  the  direction  of 
propagation . 

Finally,  we  calculate  the  mean  rate  of  energy 
dissipation  per  unit  volume  <V> .  We  could  use  the  expres¬ 
sion  for  V  in  (3.40),  but  there  is  an  easier  way.  First, 
we  write  (3.39)  in  differential  form 

(KED  +  PED)  +  V  =  -  V-I  (3.66) 

O  t 


We  have 


&  (KED+PED)  =  PUiRu.R  +  Ar0r0r  +  2MRe.jRe. 


(3.67) 


Now,  it  is  easy  to  show  that 


<Re  (g)  Re  (g)  >  =  0 


(3.68) 


where  g  =  G  (co ,  r )  elwt ,  with  G  complex 


Using  this  result. 


we  have 


. 


n 
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<  ~  (KED  +  FED) >  =  0 
0  t 


Hence,  using  (3.51), 


<V>  =  -  v  •  <I> 


->  -> 

=  2 A  •  <I> 


Evaluating  this,  using  (3.63)  to  show  that 

A  •  {(PxA)  x  (M  P  -  M  A) }  =  M  IpxAl2 

X  R  I 


(3.69) 


(3.70) 


(3.71) 


we  finally  obtain 

<V>  =  oo  |  B  |  2  e"2A’r  [  po)2  ( P  •  A)  +  4M  z  |  PxA  |  2  ]  (3.72) 

Equation  (3.70)  shows  that  the  mean  rate  of  energy  dissi¬ 
pation  per  unit  volume  is  proportional  to  the  component 
of  the  mean  energy  flux  along  A.  For  homogeneous  waves, 
we  have 


<V>  =  pco 3  |  B  |  2  e  2A*r  PA 

n 


(3.73) 


-> 


and  for  elastic  waves  (M  =0,  P*A  =  0) ,  we  have 


<V> 


0,1 


0 


(3.74) 
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3.2.2  SV  Waves 


We  now  calculate  the  same  quantities  for  SV  waves. 

t  t 
The  potential  ip  for  SV  waves  is  given  by  (3.14),  with 

k  =  P-iA  as  usual.  Using  (3.17) ,  the  displacement  is 

given  by 


u 


V  x  J 


-  -iB  m  e 


i (wt-k • r ) 


(3.75) 


where 


m  =  k  x  n 


(3.76) 


Noting  that 


0  =  V  •  u  =  V  •  (Vx-i|j)=0 


(3.77) 


equation  (3.27)  gives 

a.  .  =  -  MB  (k.m.  +  k.m.)  el(aJt"^*r)  (3.78) 

13  13  3  1 

To  obtain  the  mean  energy  flux,  we  start  by  using 
(3.41),  (3.75)  and  (3.78)  to  obtain 


<1 . >  =  -  <0 .  u  ._> 
1  13R  jR 


=  2 


2  e  2A*r  Re [ M (k.m.  +  k.m.)  m.  ] 

13  3i  3 


. 
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=  |  m  |  B  |  2  e  2A*r  Re [W { (m*m* )  k  +  (k-m*)m}]  (3.79) 


Using  (2.16),  (3.16)  and  standard  vector  identities 

involving  dot  and  cross  products,  we  get,  after  some 
tedious  algebra, 


m»m*  =  P  +  A 


'£*m*)m  =  -  2(P*A)  x  (A  +  iP) 


(3.80) 


Using  (3.80) ,  we  have 


Re  [M  (m*m* )  k]  =  [M0(P2+A2)]P  +  [Mt(PZ+A2)]a 

K  _L 


(3.81) 


From  (2.11)  we  obtain 


2  i/ 1  2 

pm  =  M  k 


(3.82) 


and  taking  the  real  and  imaginary  parts  of  this  equation, 
we  get 


d(PZ~A2)  +  2M  T ( P  *  A )  =  pm2 

R  -1 


-2M  (P  •  A)  +  M  (P2-A2)  =  0 
R  -L 


(3.83) 


We  now  rewrite  these  as 


I . 
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M_,(P2+A2)  =  pto2  +  2M  A2  -  2M,  (P  •  A) 

K  1 


(3.84) 


M  (P2+A2)  =  2M  (P • A)  +  2M  A2 
1  K  I 


Substituting  them  into  (3.81) ,  and  using  standard  vector 
identities,  we  obtain 


(3.85) 


We  also  easily  obtain 


Re[M(k*m*)m]  =  2 (P*A)  x  (M  P  -  M  A) 

X  R 


(3.86) 


Hence,  using  (3.85)  and  (3.86)  in  (3.79),  we  finally  obtain 


2A’r  [pco2  P  +  4  (PxA)  x  (M  p-MJ)]  (3.87) 

X  R 


which  has  the  same  form  as  the  mean  energy  flux  for  P 
waves.  Using  (3.70),  this  means  that  <V>  will  also  have 
the  same  form  as  the  <V>  for  P  waves . 

Only  KED  and  PED  remain.  The  mean  kinetic  energy 
density  is  easily  calculated  to  be 


<KED> 
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(3.88) 


where  (3.80)  has  been  used.  This  again  has  the  same  form 
as  the  <KED>  for  P  waves.  Hence,  we  can  conclude  that 
<PED>  will  have  the  same  form  as  <PED>  for  P  waves.  Also, 
all  of  the  results  pertaining  to  the  time-averaged  quanti¬ 
ties  for  P  waves  will  also  hold  for  SV  waves. 

3.2.3  SH  Waves 

Finally,  we  calculate  the  energy-related  quantities 


for  SH  waves .  The  potential  ip  for  SH  waves  is  given  by 


(3.22)  and  (3.23)  with  k  =  P-iA  as  always.  The  displacement 


is  given  by  (3.24)  and  (3.25).  Taking  equation  (3.77)  into 
account,  the  stress  tensor  becomes 


i  (u)t-k  •  r ) 


(3.89) 


o 


Applying  this  to  the  calculation  of  <I>,  we  get,  using 
(3.24) ,  (3.25)  and  (3.41) , 


<1 .  > 
l 


Re [M  k±] 


(3.90) 


1 


' 
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where  we  have  used  the  fact  that  k2  =  p2  “  iA2  =  °*  Hence 
<I>  =  §  OJ I  D I  2  e"2A'r  (MrP  +  MIX)  (3.91) 

We  can  rewrite  this  in  a  form  similar  to  the  P-SV  case. 
Equations  (3.81)  and  (3.85)  imply 

(P2+A2)  (MrP+MiA)  =  pco2p  +  2  (PxA)  x  (M  ?  -  M  A)  (3.92) 
hence 

i  i  2 

D  -*  0 

<I>  =  y  CO  -  e  [poi  P  +  2(PXA)  X  (M  p-M  A)  ]  (3.93) 

h  I  R 

where 


2  2 
h  =  P  +  A 


(3.94) 


Comparing  (3.93)  with  the  P-SV  result,  we  see  that  it 
differs  by  a  factor  of  2  in  the  second  vector  term.  As 
will  be  seen,  this  is  also  true  for  <PED>,  <E>  and  <V> . 

We  can  also  relate  D  to  B :  using  the  vector  identity 

|x  X  Y|  =  (XY)  -  ( X • Y )  (3.95) 

and  the  fact  that  the  condition  V  •  ip  =  0  implies 

B  •  A  =  B  •  P  =  0 


I  . 
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we  obtain,  after  some  algebra. 


+ 


-> 

B. 


2  2 
(P  +AZ) 


(B • B* )  (P2+A2) 


— > 


->  ->  * 
B  •  B 


(3.96) 


The  mean  kinetic  energy  density  can  again  be  easily  computed 
to  be 


<KED>  = 


"  2  P  UiR  UiR> 


pm 


-2A  •  r 


2  2 
(P  +  AZ) 


(3.97) 


For  the  mean  potential  energy  density,  we  again  use 
(3.77),  and  obtain 


< PED>  -  <Mr  e.  e.jR> 


1  (  [M  (P2  +  a2)2] 
4  \  h  / 


(3.98) 


We  can  again  rewrite  this  in  a  form  similar  to  the  P-SV 
case:  equation  (3.60)  implies 


Mr(P2  +  A2) 2  =  4Mr  |PxA 


+  M 


2 , 2 


R 


(3.99) 


Using 


the 


equation 


r 


we  have 
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M  =  pco^  Re  ( l/kz)  =  poo2  (p2  -  A^)/|k 


2 , 2 


R 


=o> 


MR  |  k2  |  2  =  pco 2  (P2  -  A2) 


(3.100) 


Using  this  in  (3.99)  and  (3.98),  we  finally  obtain 


<PED>  =  4- 


J 

D  j 

\  .  1 

e_2A  r  [poo2  (P2-A2)  +  4MP  |P  x  a|2] 

R 


(3.101) 

The  mean  total  energy  density  can  now  be  calculated, 
and  it  is  again  given  by  the  equation 


T  1  ->  -> 

<E>  =  -  P  •  <i> 
oo 


(3.102) 


The  mean  rate  of  energy  dissipation  per  unit 
volume  is  easily  computed  to  be 


<V>  = 


2A  •  <I> 


=  oo 


- — i  e  2A*r  [poo2  (P  •  A)  +  2  M  |  PxA  |  2] 

h  1 


(3.103) 


In  summary,  the  following  general  conclusions  can 
be  drawn : 

(a)  <I>  is  not,  in  general,  in  the  same  direction  as  P, 

except  for  elastic  waves,  and  anelastic  homogeneous 


waves 


■ 
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(b)  <PED>  ^  <KED>  except  for  homogeneous  waves  in  a 
perfectly  elastic  medium 

(c)  <E>  is  proportional  to  the  component  of  <I>  along  P 

(d)  <V>  is  proportional  to  the  component  of  <I>  along  A. 


CHAPTER  4 


REFLECTION  AND  TRANSMISSION  AT  A  BOUNDARY 

In  this  chapter,  we  treat  the  problem  of  the  reflec¬ 
tion  and  transmission  of  harmonic  plane  waves  at  a  plane 
boundary  separating  two  homogeneous  isotropic  linear 
viscoelastic  materials.  Previous  treatments  of  this 
problem  include  papers  by  Lockett  (1962) ,  Cooper  and 
Reiss  (1966),  Cooper  (1967),  Shaw  and  Bugl  (1969), 
Schoenberg  (1971)  and  Borcherdt  (1977) .  The  classical 
problem  of  a  line  source  situated  at  a  finite  distance 
from  a  plane  interface  separating  two  linear  viscoelastic 
solids  has  been  examined  by  Buchen  (1971b) . 

The  treatment  of  plane  wave  reflection  and  trans¬ 
mission  by  Borcherdt  (1977)  is  the  most  recent  and  the 
most  general,  and  provides  part  of  the  theoretical  frame¬ 
work  for  the  computation  of  reflection  and  transmission 
coefficients  in  this  chapter.  In  the  first  section,  we 
treat  the  problem  of  the  reflection  and  transmission  of  SH 
waves  in  detail,  and  examine  some  of  the  interesting 
characteristics  exhibited  by  viscoelastic  waves  which  are 
not  exhibited  in  the  perfectly  elastic  case.  We  then 
discuss  the  computation  of  reflection  and  transmission 
coefficients  and  present  some  examples  to  illustrate  the 
differences  between  elasticity  and  viscoelasticity.  The 
treatment  of  the  P-SV  case  is  very  similar  to  the  SH  case, 
and  is  outlined  in  the  second  section. 
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4 . 1  SH  Waves 

Figure  1  shows  the  incident,  reflected  and  trans¬ 
mitted  waves  at  a  boundary  separating  two  viscoelastic 

t  — y 

media  V  and  V ' .  The  complex  propagation  vectors  k  for  the 
incident  and  reflected  waves  can  be  written  as 

/\  •  /\ 

kgj  =  kSxx  +  (-l)3dg  z  (3=1,2)  (4.1) 


where 


d 


3 


±p.  v. 


(4.2) 


and  where  j=l  denotes  the  incident  wave,  and  j=2  the 

/\  /\ 

reflected  wave,  and  x  and  z  are  unit  vectors  in  the  x  and 
z  directions.  For  the  transmitted  wave,  we  have 


(4.3) 


where 


±p.  v 


(4.4) 


The  symbol  p.v.  denotes  the  principal  value  of  the  complex 
square  root,  i.e.,  the  complex  root  for  which  -90°<ri£  +  90  , 
where  n  is  the  angle  that  the  root  makes  in  the  complex 
plane.  The  choice  of  sign  for  d^  and  d^  will  be  discussed 
further  on  in  this  section. 

In  these  equations,  we  have  set 

kSx(inc.)  =  k  Sx(ref  1 .  )  =  kSx(transm.)  =  kSx 


(4.5) 


. 
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FIGURE  1:  Incident,  reflected  and  transmitted  rays  at  a 
boundary  separating  two  anelastic  media.  The 
angles  y  lie  between  -90°  and  +  90°.  They  are 
positive  as  shown. 
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in  anticipation,  since  this  will  be  the  condition  required 
to  satisfy  the  boundary  conditions  of  this  problem. 

From  the  equations  kg  =  Pg  -  iA  ,  we  can  calculate 
the  progation  and  attenuation  vectors  for  the  three  waves. 
Suppressing  the  subscript  "S"  from  hereon  in,  we  obtain 


P  . 
3 


A  . 
3 


kxR  X  +  (-1):l  dBR  Z 
-kxi  X  +  (-D3+1a6I  Z 


j  =  l ,  2 


(4.6) 


P  ' 


=  k 


xR 


x  - 


d ' 

SR 


A'  =  -k  _  x  +  d '  z 
xl  SI 

where  the  subscripts  "R"  and  "I"  stand  for  the  real  and 
imaginary  parts,  as  before. 

To  obtain  the  reflection  and  transmission  coeffi¬ 
cients  (by  which  we  mean  the  complex  relative  wave 
amplitudes  of  reflection  and  transmission,  and  not  the 
squares  of  their  magnitudes )  ,  we  assume  the  two  media 
are  in  welded  contact,  which  means  we  require  the  stress 
and  displacement  to  be  continuous  across  the  boundary. 
The  stress  and  displacement  for  SH  waves  are  given  by 
(3.89)  and  (3.24)  respectively.  Hence,  we  have  at  z=0. 


u  =  u' 


a 


32 


33 


=  0) 


(4.7) 


where 
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u 


O  ■  /  7*"  . 

2  i  (cot-k  .  •  r )  ^ 

Z  D  .  e  3  y 

j  =  l  : 


(4.8) 


i  (cot-k '  *  r ) 


u'  =  D'  e 


y 


Assuming,  for  the  moment,  that  the  three  k  ' s  are  unequal, 
equations  (4.7)  lead  to 


-ik  ..x  -ik  „x  -ik'x 

Dxe  xl  +  D2e  x2  =  D'e  x 

“ikxlx  _ikv?x  -ik'x 

Md6iDie  -Md02D2e  =  M'd'D'e  x 


(4.9) 


Since  the  complex  amplitudes  D^,  and  D'  are  independent 
of  the  spatial  coordinates  x,  v  and  z,  equation  (4.9) 
implies 


k  =  k  =  k '  =  k 
xi  x2  x  x 


(4.10) 


which  is  the  equivalent  of  equation  (4.5).  Equation  (4.5) 
(or  (4.10))  leads  to  the  anelastic  form  of  Snell's  lav/,  as 
we  shall  see  below.  Equations  (4.9)  now  become 


D1  +  °2  ’  D' 


Mdg(D1-D2)  =  M'd^D’ 


and  the  solutions  are 


(4.11) 


D. 


D. 


D' 

D, 


Md3  -  M'd' 
Mdg  +  M'd£ 


(4.12) 


2Md 


6 


Md^  +  M'd£ 


.  i 


■ 


The  first  equations  in  (4.12)  gives  the  reflection 
coefficient,  and  the  second  gives  the  transmission 
coefficient . 


The  SH  reflection  coefficient  at  a  free  surface  is 

easily  obtained  from  (4.11)  by  equating  the  stress  at  the 

surface  to  zero,  which  yields  D]_=D2'  i*e*'  (D2/D^)  =  ^ * 

As  mentioned  above,  equation  (4.5)  leads  to  Snell's 

law  for  viscoelastic  media.  Since  k  =  P  -iA  ,  there  are 

x  xx 

now  two  parts  to  Snell's  law,  i.e.  A  ,  as  well  as  P  ,  must 

x  x 

be  continuous  across  the  boundary.  Referring  to  Figure  1, 
this  implies 


P,  sin©..  =  P„  sin0o  =  P'  sin0 '  =  k  _ 
1122  xR 

A1  sin(01-Y  )  =  A2  sin ( 0 2“Y2 )  =  A'  sin(0'-y') 


(4.13) 


This  is  the  anelastic  version  of  Snell's  law.  If 

■  ->  -  . 

|  v  I  =  a)/ 1 P  |  is  substituted  into  the  first  equation  of  (4.13) 
it  takes  the  form 


sin0^ 

sin0  2 

sin0  ' 

i  i 

1  1 

i  -*■ ,  | 

1  V!  1 

lv2l 

1  v '  | 

(4.14) 


where  v. 


V. 


and  v'  are  the  phase  speeds  of  the 


1 1  '  1  v2 

incident,  reflected  and  transmitted  waves.  From  Chapter  2, 
we  see  that  the  phase  speeds  depend  on  y^,  y2  and  y' 
respectively,  and  since  y'  will  generally  depend  on  the 
angle  of  incidence  0^  the  phase  velocity  of  the  trans- 
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mitted  wave  will  also  depend  on  0  . 

Equations  (4.6)  show  that  |p  |  =  |P  |  and 
| Ail  =  | A2 |  ,  hence  equations  (4.13)  show  that  0 ^=9 2  and 
Yl=Y2- 

A  critical  angle  for  the  transmitted  SH  wave  is 
defined  as  an  angle  of  incidence  for  which  the  transmitted 
SH  wave  propagates  parallel  to  the  interface.  The 
following  theorems,  originally  due  to  Borcherdt  (1977) , 
show  that  critical  angles  for  viscoelastic  media  differ 
in  nature  from  those  for  elastic  media. 

Theorem  1:  If  medium  V  is  perfectly  elastic,  and  medium  V' 
anelastic,  then  there  are  no  critical  angles  for  the 
transmitted  SH  wave. 

Proof:  Suppose  0^  is  a  critical  angle  for  the  transmitted 
SH  wave,  i.e.  0  '  =tt/2  ,  and  P'  is  in  the  x  direction.  Since 
V  is  perfectly  elastic,  A^=0  and  (4.13)  shows  that  k^I=0. 
Hence,  from  (4.6),  A'  is  in  the  z  direction,  which  means 
y'  =  ±  tt/2,  i.e.  P'*A'=0.  But  we  know  from  Chapter  2  that 
such  a  wave  cannot  propagate  in  an  anelastic  medium.  Hence, 
there  are  no  critical  angles  in  this  case. 

Theorem  2:  If  medium  V  is  anelastic  and  0^  is  a  critical 
angle  for  the  transmitted  SH  wave  (0^tt/2),  then 

sin  01  sin(01~y1)  =  £  cosy1  (4.15) 


where 
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r  <k'2)I  kRkl 

‘  (k* 2)I  _  kRkI 

Proof :  If  0^  is  a  critical  angle 
then  the  equation  for  p'  in  (4.6) 
implies 


(4.16) 


for  the  transmitted  SH  wave, 
implies  d'  =  0,  which 

P  K 


2deRdei  “  Im[dB2]  =  (k'2)i  -  <k2)i  “  0 


(4.17) 


Now  using  (4.13)  and  (2.16),  we  have 


(kx)x  lm[ (kxR  +  lkxi)  ^  2kxRkxI 


=  -2P^A^  sin  0^  sin(0  -y  ) 

(k2)x 

= -  sin  0.  sin(01-y_) 

cosy.  1  1  '1 


Using  this  in  (4.17)  gives  the  required  relation  (4.15). 
The  above  theorem  is  stated  and  proved  in  a  form  slightly 
different  from  that  of  Borcherdt  (1977) . 


Corollary  2.1:  If  V  is  anelastic  and  V'  elastic,  then  there 

exists  at  most  one  critical  angle,  namely  0^=y^. 

2 

Proof :  If  V'  is  elastic,  then  ( k *  )  =  0,  hence  £>0,  and 

(4.15)  becomes 

sin  0^  sin(0^-y^)  =  0 

Hence,  if  0^  is  critical,  then  0^=y^  . 


(4.18) 


3 


■ 
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Theorem  2  shows  that  critical  angles  can  exist  only 

when  (4.15)  is  satisfied,  i.e.,  for  distinct  values  of 

0^  and  This  means  that  critical  angles  in  viscoelastic 

media  are  discrete.  In  other  words,  if  0  is  a  critical 

c 

angle  and  0.  is  increased  beyond  0  ,  then  0  is  no  longer  a 

C  _L 

critical  angle,  i.e.  the  transmitted  propagation  vector  p' 
is  no  longer  parallel  to  the  interface.  p'  is  parallel  to 
the  interface  only  for  0^=0^.  This  result  represents  a 
significant  deviation  from  perfect  elasticity  theory. 

Cooper  (1967)  obtained  the  same  result.  Note  that  Theorem  2 
does  not  imply  that  every  solution  0^  of  (4.15)  is  a 
critical  angle;  0^  could  be  an  angle  for  which  A',  rather 
than  ,  is  parallel  to  the  interface. 

For  two  anelastic  media,  the  quantities  pertaining 
to  the  transmission  medium,  i.e.  P',  A',  0'  and  y',  can  be 
calculated  as  follows.  First  of  all,  we  assume  that  0^ 
and  y^  are  initially  known,  along  with  the  medium  parameters, 
i.e.  the  densities,  homogeneous  phase  speeds,  and  loss 

2 

factors.  P,  ,  and  hence  k  can  then  be  calculated.  k 
11  x 

2 

and  k'  can  be  calculated  solely  from  the  medium  parameters. 
We  can  then  see  from  (4.6)  that  the  magnitudes  of  the 
transmitted  propagation  and  attenuation  vectors  P'  and  A' 
are  given  by 


P ' 


A' 


+  d 


+  d 


(3R 


31 


(4.19) 
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The  angle  of  transmission,  0* ,  is  obviously  given  by 

"1  ^xR 

0'  =  tan  (-3^)  (4.20) 

3R 

Once  P',  A'  and  0'  are  known,  y1  can  then  be  calculated 
from  (4.13).  Since  Q  is,  in  general,  frequency-dependent, 
velocities  and  angles  will  also,  in  general,  depend  on 
frequency . 

Equation  (4.20)  shows  that  0'  depends  on  d'  ,  and 

pR 

as  mentioned  earlier,  there  is  a  choice  of  sign  to  be 
considered  in  the  calculation  of  d^,  i.e.,  a  choice  of 
which  of  the  two  square  roots  to  use  (see  (4.4)).  Intuitive¬ 
ly,  one  would  expect  that  the  "+"  sign  would  always  be 

chosen,  because  this  would  make  d'  always  positive,  and 

pR 

hence  from  (4.6)  P'  would  always  point  in  the  proper 
direction,  i.e.  it  would  always  point  into  the  transmission 
medium.  However,  such  a  choice  of  sign  has  been  found 
to  introduce  errors  in  the  computation  of  the  reflection 
and  transmission  coefficients  at  and  beyond  the  critical 
angle  (if  one  exists).  For  example,  assuming  that  a  critical 
angle  ©c  for  the  transmitted  SH  wave  exists,  let  us 
consider  an  angle  of  incidence  0^  slightly  less  than  0^, 
and  suppose  that  d^,  for  this  angle  of  incidence  0^,  lies 
in  the  4th  quadrant  of  the  complex  plane,  near  the  imaginary 
axis.  In  other  words,  d'  makes  an  angle  slightly  greater 

p 

than  -90°  in  the  complex  plane,  i.e.  d£  is  almost  purely 
imaginary  (it  is  often  the  case  that  d^p^O  and  d^^O  for 


' 


s 


■ 


P'  is  in  the  x 


0^0, <9  ) .  For  the  critical  angle  0  , 

1  c  c 

direction,  i.e.  parallel  to  the  interface,  hence  (4.6) 

implies  dgR=0  and  that  dg  is  purely  imaginary,  i.e., 

/  2  2^ 

dg  =  +  p. v. /  k*  -k^  =  +  ia,  where  "a"  is  a  constant.  If 
we  choose  the  "+"  sign,  i.e.  the  "positive"  root,  then, 
due  to  the  definition  of  the  principal  value,  d^  undergoes 
a  sudden  change  of  sign,  from  negative  for  0  c 0 c  r  to 
positive  for  0=0  .  As  we  can  see  from  (4. 12), this  would 

JL  C 

correspondingly  introduce  a  sudden  change  in  the  reflection 

and  transmission  coefficients  at  the  critical  angle,  i.e. 

a  jump  discontinuity  at  0  (see  Figure  2a)  To  avoid  this, 

o 

we  must  choose  the  sign  at  0^.  For  super-critical 

2 

incident  angles,  the  same  problem  occurs  since  dg  ,  which 
lies  along  the  negative  real  axis  of  the  complex  plane 
for  0=0  ,  moves  upwards,  in  this  particular  example,  into 

-L  C 

the  2nd  quadrant  for  0,>0  ,  which  means  that  choosing  the 

"+"  sign  for  dg,  i.e.  the  "positive"  root,  would  again 

cause  the  same  sudden  change  of  sign  (-  to  +)  in  dg^,  due 

to  the  definition  of  the  principal  value.  Hence,  to 

eliminate  the  jump  discontinuities  in  the  reflection  and 

transmission  coefficients,  i.e.  to  obtain  smooth  continuous 

curves,  we  must  choose  the  sign  of  dg  appropriately.  In 

this  example,  choosing  the  "  +  "  sign  for  0,<0c  an(i  the 

sign  for  0,^0  will  ensure  that  d'  will  pass  smoothly  from 
1  c  p 

the  4th  quadrant  for  01<0C  into  the  3rd  quadrant  for  01*0c 
(rather  than  the  1st) ,  of  the  complex  plane,  and  hence  no 


. 


■ 
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FIGURE  2:  SH  anelastic  reflection  coefficients  illustra¬ 
ting  jump  discontinuities.  The  top  curve  in 
each  diagram  is  the  relative  amplitude  and  the 
bottom  is  the  relative  phase.  The  loss  factors 
for  the  media  of  incidence  and  transmission  are 
RQ1  and  RQ2 ,  respectively.  DENI  and  DEN2  are 
the  densities  in  g/cc .  VH1  and  VH2  are  the 
homogeneous  phase  speeds  in  km/s .  The  data, 
taken  from  Silva  (1976),  represent  highly  atten- 
uative  soil  layers.  The  incident  attenuation 
angle  y^  is  40°  for  all  three  diagrams.  The 
jump  discontinuities  occur  because  of  a  wrong 
choice  of  sign  for  d^': 

(a)  the  positive  principal  value  in  equation 
(4.4)  is  chosen  for  all  angles  of  incid¬ 
ence.  A  jump  discontinuity  is  produced  at 
the  critical  angle  0^  =  55.53°. 

(b)  The  sign  of  d^'  is  chosen  by  requiring  the 
mean  energy  flux  of  the  transmitted  ray  to 
always  point  away  from  the  boundary,  for  all 
angles  of  incidence.  The  jump  discontinuity 
occurs  at  the  angle  of  incidence  for  which 
the  mean  energy  flux  of  the  transmitted  ray 
is  parallel  to  the  interface. 

(c)  The  sign  of  d^'  is  chosen  correctly,  hence 
no  jump  discontinuity  is  produced  and  the 
curves  are  continuous . 
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jump  discontinuities  will  occur  (see  Figure  2c) . 

The  same  arguments  can  be  applied  when  d^^O  for 
O£0^<0c.  In  these  cases,  we  must  choose  the  "+"  sign  for 
0, £0  and  sign  for  0_>0  . 

Alternatively,  Borcherdt  (1977)  has  suggested  that 
the  signs  of  d^  and  d^  be  determined  by  requiring  that  the 
mean  energy  fluxes  associated  with  the  incident  and  trans¬ 
mitted  waves  point,  respectively,  toward  and  away  from  the 
boundary.  From  a  theoretical  viewpoint,  this  seems  like  a 
very  natural,  physically  realistic  criterion.  However,  as 
it  turns  out,  this  criterion  also  produces  jump  discontinui¬ 
ties  in  the  reflection  and  transmission  curves,  although 

— y 

not  necessarily  at  the  critical  angle,  for  which  P'  is 

parallel  to  the  interface,  rather,  at  the  angle  at  which 

— 

the  mean  energy  flux  of  the  transmitted  wave,  <I'>,  is 

parallel  to  the  interface  (see  Figure  2b) . 

Apparently,  we  must  accept  the  criterion  proposed 

above,  involving  the  reversal  of  the  sign  of  d^  at  0^=0c. 

However,  while  this  criterion  does  eliminate  the  jump 

discontinuities,  is  also  introduces  conceptual  problems 

involving  the  transmitted  wave.  For  instance,  the  angle  of 

transmission  is  given  by  (4.20).  However,  when  0-^>0c,  we 

have  d'  <0,  according  to  the  above  criterion,  hence  0'  is 
pR 

apparently  negative.  Actually  this  is  not  the  case:  what 
we  must  do  is  use  the  branch  of  the  arctangent  function 
immediately  above  the  principal  branch,  i.e.. 


. 
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6'  =  tan  X(-^)  =  tt  +  p.v.  [tan  1  (-^)  ]  (4.21) 

pR  3R 

where  p.v.  is  the  principal  value  of  the  arctangent,  i.e. 

the  value  on  the  principal  branch.  This  means  that  for 

9-^>0c^  we  have  0'>9O°.  This  interpretation  is  consistent 

with  equation  (4.6)  which  shows  that  P'  points  into  the 

incidence  medium  when  d'  <0.  This  seemingly  strange 

situation  must  be  accepted  if  we  are  to  avoid  jump 

discontinuities  at  0  in  the  reflection  and  transmission 

c 

coefficients.  We  have  already  shown  in  the  discussion  of 
Theorem  2  that  P'  must  move  away  from  the  interface  for 
0^0  ,  and  so  the  situation  is  not  really  that  surprising. 

The  only  surprise  is  that  P'  must  move  into  the  incidence 
medium  rather  than  back  into  the  transmission  medium. 

However,  it  usually  does  not  move  very  far  into  the  incidence 
medium,  i.e.  0'  usually  exceeds  90°  by  only  a  very  small 
amount,  on  the  order  of  a  fraction  of  a  degree,  on  the 
average . 

One  might  think  that  this  situation  occurs  because 
of  the  plane-wave  assumption.  However,  Buchen  (1971b)  has 
reached  the  same  conclusions  in  his  treatment  of  the 
classical  problem  of  a  line  source  radiating  cylindrical 
waves,  situated  at  a  finite  distance  from  a  plane  inter¬ 
face  separating  two  viscoelastic  media.  He  suggests  that 
no  true  ray  path  can  be  associated  with  this  transmitted 
wave  beyond  the  critical  angle  because  it  does  not  satisfy 
Fermat's  principle. 


- 
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One  should  interpret  the  transmitted  wave  with 
Q'^O0  not  as  an  actual  ray,  but  rather  as  representing 
a  small  energy  flux  across  the  interface  from  V'  to  V 
for  6]_>9C*  This  does  not  occur  in  the  perfectly  elastic 
case.  Borcherdt  (1977)  has  obtained  results  concerning 
energy  flow  across  the  interface  which  support  this  notion. 
He  has  shown,  for  SH  viscoelastic  waves,  that  the  normal 
component  of  the  total  mean  energy  flux  (i.e.,  the  flux  for 
the  total  wave  field)  is  continuous  across  the  plane 
boundary  separating  V  and  V' .  The  proof  of  this  rests  on 
the  fact  that  the  z  component  of  the  mean  energy  flux  can 
be  shown  to  be  proportional  to  cr^'  which,  according  to 
(4.7),  must  be  continuous  across  the  boundary.  In 
mathematical  terms,  we  have 


«I  >  +  <I2>  +  <^12>  +  <*21>^*z  =  <*'>*z 


(4.22) 


where  <I^>  and  <l2>  are  t^ie  mean  energy  fluxes  associated 


with  the  incident  and  reflected  waves,  <I’>  is  the  mean 
energy  flux  associated  with  the  transmitted  wave,  and 
<1  2>  ( < 1 2i> )  is  the  mean  energy  flux  due  to  the  velocity 

field  of  the  incident  (reflected)  wave  interacting  with  the 
stress  field  of  the  reflected  (incident)  wave.  <I12>  and 
<1  >  do  not  occur  in  the  perfectly  elastic  case.  If  we 


define 


•  ->  .  ~y 

RE  |  <  X  2  >  *  Z/<I]_>  *  z 


T  E  <1 ’ > • z/<I1> • Z 


IC. .=  -  <1. .>«z/<J- >*z 

ID  ID  1 


(4.23) 
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then  equation  (4.22)  can  be  written  as 


R  +  T  +  IC12  +  IC21  =  1 


(4.24) 


which  states  that  energy  is  conserved  at  the  boundary.  In 
the  perfectly  elastic  case,  IC^2  +  IC21  =  ^  anc^  (4.24) 
simplifies  to  the  familiar  form  R  +  T  =  1. 

Theorem  2  shows  that  it  is  possible  that  two  critical 
angles  for  the  transmitted  SH  wave  may  exist,  since  (4.15) 
can  have  two  solutions  for  0^.  Numerical  computations 
verify  that  such  situations  do,  in  fact,  exist  in  certain 
instances.  In  these  cases,  the  transmitted  angle  0'  behaves 
in  the  following  way  for  monotonically  increasing  0^:  let 


0C^  and  0^o  be  the  two  critical  angles,  such  that  0r1<0pO. 


c2 


cl  c2 


(o 


For  0^<0C-^,  we  have  0'<9O  .  After  the  first  critical  angle 

0  ,  is  achieved,  0'  exceeds  90°,  reaches  a  maximum  value, 

^  1 

decreases  back  to  90°  for  the  second  critical  angle  0 c2, 

and  then  continues  to  slowly  decrease,  remaining  near  90°, 

until  0^  is  parallel  to  the  interface. 

Concerning  the  choice  of  sign  for  d^  in  equation 

(4.2),  we  assume  that  0.^  always  satisfies  0^0^90°.  Hence 

(4.6)  implies  d^R  is  always  >0.  This  means  that  the  "+" 

sign  is  always  chosen  for  d^. 

We  can  see  from  the  above  discussions  that  if  we 

restrict  our  calculations  to  sub-critical  regions,  i.e. 

0  <0  ,  then  the  calculations  are  substantially  simplified, 

1  c 

since  we  do  not  need  to  worry  about  choosing  the  correct 
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sign  for  d^,  or  computing  0'  correctly,  etc.  However,  to 
obtain  a  more  thorough  understanding  of  viscoelastic 
reflection  and  transmission,  it  is  necessary  to  examine  the 
super-critical  behaviour  of  viscoelastic  waves. 

In  the  following  chapters,  synthetic  seismograms  will 
be  computed  for  layered  anelastic  media.  A  ray  theory 
approach  will  be  used,  which  means  that,  for  a  given  ray, 
it  will  be  necessary  to  compute  reflection  and  transmission 
coefficients  at  all  the  interfaces  which  the  ray  encounters. 
In  order  to  be  able  to  compute  synthetic  seismograms  for 
a  general  layered  medium,  which  may  contain  perfectly 
elastic  layers  embedded  among  anelastic  layers,  we  must  be 
able  to  compute  reflection  and  transmission  coefficients 
for  two  half-spaces  separated  by  a  boundary,  where  each 
half-space  may  be  either  elastic  or  anelastic.  In  order  to 
do  this  effectively,  we  divide  the  calculations  of  reflec¬ 
tion  and  transmission  coefficients  into  six  separate  cases: 

anelastic  incident  rays 
anelastic  transmitted  rays 

anelastic  incident  rays 

inhomogeneous  elastic  ("JEL") 
transmitted  rays  (y'=+90  ) 

IEL  incident  rays  (y^=+90°) 
anelastic  transmitted  rays 


(1)  fQ  Vo  , 

\q'~Vo  , 

(2)  fcfVo  , 

Q' _1=0  , 
Q~1=0  , 

1q,‘Vo  , 


(3) 
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|  -i 

(4) Q  =0  ,  IEL  incident  rays 

'  ‘*’=0  ,  IEL  transmitted  rays 


(5)  |Q  1=0 


elastic  incident  rays 
anelastic  transmitted  rays 


(6)  Q  ■*■=()  ,  elastic  incident  rays 

\  -i 

Q'  =0  ,  elastic  transmitted  rays 


Case  1  involves  two  anelastic  media,  and  critical 


angles  can  be  determined  by  solving  (4.15)  for  0^  and 
checking  to  see  if  0  is  a  critical  angle  (i.e.  if  6'=90°). 


Case  2  shows  that  if  V  is  anelastic  and  V'  elastic 


the  transmitted  rays  must  be  inhomogeneous  elastic  ("IEL") 
rays  (see  Figure  3).  This  can  be  demonstrated  as  follows. 
Since  A^  sin(0  -y^)  is  not,  in  general,  equal  to  zero, 
equation  (4.13)  implies  that  A*  is  also  not  zero,  in  general, 
and  since  V'  is  elastic,  we  know  that  P'  and  A'  must  be 
perpendicular  to  each  other.  Also,  A  must  be  continuous 

X 

across  the  interface.  Hence,  the  general  situation  is  shown 
in  Figure  3.  We  see  that  the  continuity  of  A^  can  lead  to 
a  situation  in  which  the  transmitted  wave  is  damped  toward 
the  interface,  rather  than  away  from  it.  This  can  also 
happen  in  the  other  cases.  Such  situations  have  also  been 
found  by  other  authors  such  as  Cooper  and  Reiss  (1966), 

Cooper  (1967)  and  Buchen  (1971b).  Corollary  2.1  shows 
that  if  a  critical  angle  exists,  it  satisfies 


■ 
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FIGURE  3:  Ray  diagram  illustrating  the  appearance  of  an 
inhomogeneous  elastic  ( " IEL" )  transmitted  ray 
in  V',  for  the  case  in  which  V  is  anelastic  and 


V'  elastic. 


4 
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Cases  3  and  4  are  treated  in  order  to  cover  the 
possibility  of  the  transmitted  IEL  ray  of  Case  2  in  turn 
being  an  incident  ray  for  another  interface.  These  cases 
were  not  treated  by  Borcherdt  (1977) . 

For  case  3,  we  can  prove  a  theorem  concerning 
critical  angles  similar  to  Theorems  1  and  2. 

Theorem  3;  Suppose  V  is  an  elastic  medium  in  which  IEL 
incident  rays  (y^=+90°)  propagate,  and  V'  is  anelastic. 

Then,  if  y^=+90°,  there  are  no  critical  angles.  If  y^=-90°, 
and  0^  is  a  critical  angle  for  the  transmitted  SH  wave,  then 
6^  satisfies 

2 

(k '  ) 

sin  20t  =  -  - —  (4.25) 


Proof:  Assume  0^  is  a  critical  angle  for  the  transmitted 
wave.  Then  (4.6)  implies  d^R=0,  which  implies 


2d6Rdh  =  Im(d62)  -  (k'2)I  *  (kx>I  =  0 


(4.26) 


,o 


Now,  using  (4.13),  (2.16)  and  the  fact  that  y1=+90  ,  we 

have 


(k  ) _  =  2k  k 
x  I  xR  xl 


=  “2P^A^  sin0^  sin[0^  -  (+90°)] 

=  +  2P^A^  sin©1  cos01 

=  +  P1A1  sin  20  (4.27) 


Hence,  (4.26)  implies 


. 


must  lie  in 
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sin20_,  =  + 


(k'* 2) 


1  P1A1 

If  y1=+90°,  then  since  (k'2)i  is  negative,  0  must  lie  in 
the  range  7r/2£0^£7T,  which  is  unphvsical.  So  0^  cannot  be  a 
critical  angle  if  y1=+90°.  For  y  =-90°,  we  obtain 
equation  (4.25) . 

Note  that  Theorem  3  shows  that  there  can  possibly  be 
two  critical  angles  for  the  transmitted  SH  wave,  since 
(4.25)  yields  two  solutions  for  0  ,  as  for  Theorem  2. 

For  Case  4,  the  necessary  continuity  of  A  across 
the  interface  shows  that  incident  rays  which  are  IEL,  in 
general  must  produce  transmitted  rays  which  are  also  IEL. 
Concerning  critical  angles  for  Case  4,  we  have 
Theorem  4:  If  V  and  V'  are  elastic  media  in  which  IEL  rays 

propagate,  then  there  are  no  critical  angles  for  the 
transmitted  SH  wave. 

Proof :  Assume  0^  is  a  critical  angle  for  the  transmitted 
SH  wave.  Hence  (4.6)  implies  d'  =0,  which  in  turn  implies 

pK 

2 

equation  (4.26).  Since  V’  is  elastic,  we  have  (k'^)  =0. 

2 

Equation  (4.26)  thus  implies  (k  )  =0.  Hence  (4.27)  then 

X  1 

yields  9^=0,  tt/2  ,  ...,  which  is  unphysical  for  a  critical 

angle.  Therefore,  no  critical  angles  exist. 

For  Case  5,  Theorem  1  shows  that  there  are  no 
critical  angles  for  the  transmitted  SH  wave. 

Case  6  is  simply  the  case  for  two  perfectly  elastic 


half-spaces . 


i  p 


KH 
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Cases  1-6  are  all  the  possible  cases.  For  instance, 
the  case 


Q  =0  ,  IEL  incident  rays 

Q'  ^=0  ,  elastic  transmitted  rays 


is  NOT  possible.  The  transmitted  rays  have  to  be  IEL  rays. 

In  order  to  get  an  idea  of  what  anelastic  reflection 
and  transmission  coefficients  look  like,  we  plot  some 
examples  of  anelastic  reflection  and  transmission  amplitudes 
and  phases.  Figures  4-7  are  plots  of  reflection  and 
transmission  coefficients,  the  top  graph  being  the  amplitude, 
and  the  bottom  graph  being  the  phase,  for  SH  waves  impinging 
upon  a  boundary  separating  two  anelastic  media  1  and  2. 

The  parameters  for  each  medium  are  the  S  wave  loss  factor 
Q  \  the  density  P,  and  the  homogeneous  phase  speed  vu(  and 
they  are  shown  at  the  top  of  each  figure.  The  densities 
(DENI  and  DEN2)  are  in  g/cc  and  the  homogeneous  phase  speeds 
(VHl  and  VH2)  are  in  km/s .  The  medium  parameters  used  are 
the  parameters  of  the  first  two  layers  of  a  crustal  model 
used  to  compute  synthetic  seismograms  in  the  following 
chapters.  The  coefficient  notation  is  such  that  H1H2,  for 
example,  denotes  the  transmission  coefficient  for  an  incident 
SH  ray  in  medium  1  and  a  transmitted  SH  ray  in  medium  2. 

Each  plot  shows  curves  (#2,  3  and  4)  for  different 

values  of  Y .  ,  as  well  as  the  curve  (#1)  for  the  perfect- 

ly  elastic  case  (y.  denotes  the  value  of  the  attenuation 
1  me 

angle  for  the  incident  ray  -  it  corresponds  to  y^  in 


FIGURE  4: 


HlHl  reflection  coefficients.  Curve  #1  is  for 


the  perfectly  elastic  case.  Curves  #2,  3  and 
4  correspond  to 

(a)  Yj_nc  =  -40°,  -60°  and  -80°,  respectively 

(b)  Y^nc  =  20° ,  0°  and  -20°,  respectively 

(c)  Yj_nc  =  40°,  60°  and  80°,  respectively 

where  y.  is  the  attenuation  angle  for  the 
incident  ray.  The  phase  curves  correspond  to 
amplitude  curves  #3,  i.e.,  to  y .  =  -60°,  0° 

and  60°  in  diagrams  (a),  (b)  and  (c)  respec¬ 

tively  . 
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H1H1  RQ 1  =  0.03333  RQ2=  0.02222 

DENI =  2.100  DEN2=  2.600 

VH1=  2 .400  VH2=  3.500 


ANGLE  OF  INCIDENCE  IN  DEGREES 


0  1  0  20  30  40  50  60  70  80  90 

ANGLE  OE  INCIDENCE  IN  DEGREES 


-t— I — I — lilt) — — I  I  I  I  I  I  I  HH — I — III — III  — I  I  I  I — I — I — I-  +■  I-  f-+-|  I — I- 


<XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx 


,X* 


XXX 


xxxx 


xxxxxx 


XXXXXXXXXXXXXXXXX 


xxxxxxxxxw 


X  vx' 
xxxXX 


I  I  I  I  I  I  I  I  I  I  I  I  I  H  I-  I  I  I  I  I  I  I  I  I  I  I  1  I  I  I  I  1  I  I  I-  I  t  I  I  I  I 


(c) 


FIGURE  5 


H1H2  transmission  coefficients, 
see  caption  of  Figure  4 . 


For  details, 
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FIGURE  6 


H2H2  reflection  coefficients, 
see  caption  of  Figure  4 . 
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FIGURE  7 


H2H1  transmission  coefficients.  For  details, 
see  caption  of  Figure  4 . 
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Figure  1) .  In  each  plot,  the  phase  curve  corresponds  to 

amplitude  curve  #3.  Only  one  phase  curve  has  been  plotted 

in  each  diagram  to  avoid  unnecessary  cluttering.  The  three 

phase  curves  for  each  reflection/transmission  coefficient 

correspond  to  y.  =-60°,  0°  and  +60°  and  illustrate  the 

me 

general  behaviour  of  the  phases  as  a  function  of  Y . 

^  me 

When  the  incident  ray  is  in  medium  1  (HlHl  and  H1H2, 

i.e.  Figures  4  and  5)  critical  angles  for  the  transmitted 

ray  occur  only  for  values  of  Y^nc  greater  than  or  equal 

to  about  +17.43°.  For  lesser  values  of  y,  ,  the  trans- 

'  me 

t  t  —y 

mitted  propagation  vector  (P  *  in  Figure  1)  never  quite 

becomes  parallel  to  the  interface.  No  critical  angles 

occur  when  the  incident  ray  is  in  medium  2.  As  we  have 

previously  seen,  the  numerical  values  of  these  critical 

angles  vary  with  Yj_nc*  For  the  perfectly  elastic  case, 

the  critical  angle  for  the  transmitted  ray  is  43.29°. 

For  the  anelastic  case,  critical  angles  occur  in  Figures 

4b,  4c,  5b  and  5c  for  y.  =20°,  40°,  60°  and  80°,  and  are 
'  '  'me 

44.74°,  56.70°,  69.63°  and  83.14°  respectively.  At 

y.  =y  &17.43°,  the  value  of  y.  for  which  critical  angles 
'me  'nr'  me 

begin  to  appear,  the  critical  angle  has  the  same  value  as 

for  the  perfectly  elastic  case,  i.e.  0  (y  )=0  (elastic)= 

^  -1  c  m  c 

43.29°.  The  exact  value  of  y  can  be  determined  by 

m 

substituting  sin  0 . =v„/v '  (the  formula  for  the  critical 

_L  H  H 

angle  in  the  perfectly  elastic  case)  into  (4.15).  One 


obtains 


. 


V.*’ 
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tan  y  =  ( 1-G) tan©., 
m  1 


%/v') (1-G) 


(4.28) 


/l- (VR/V') 2 


where 


G  =  (v£/v  ) /  =  ^ 


Q 


(4.29) 


With  the  medium  parameters  shown  in  Figures  4-7,  i.e. 

vTT  =  2.4  km/s,  v'=3.5  km/s ,  etc.,  we  get  y  =17.42866...° 
H  H  m 


m 


We  can  see  that  unlike  the  perfectly  elastic  amplitude 
curve  which  exhibits  a  sharp  transition  when  going  from 
sub-critical  to  super-critical  incident  angles,  the 
anelastic  amplitude  curves  show  a  smooth,  rapid  transition. 
This,  combined  with  the  discrete  nature  of  critical 
angles  shown  previously,  shows  that  in  linear  viscoelastici¬ 
ty,  the  concept  of  a  critical  angle  is  somewhat  vague  or 
degraded.  Other  authors,  for  example,  Cooper  (1967), 

Buchen  (1971b),  Kennett  (1975),  Borcherdt  (1977)  and 
Fryer  (1978)  have  arrived  at  similar  conclusions.  In 
Figures  6  and  7,  where  the  incident  ray  is  in  medium  2  (no 
critical  angles  occur)  the  curves  do  not  deviate  from  the 
perfectly  elastic  curve  (#1)  as  much  as  in  Figures  4  and  5. 

4.2  P-SV  Waves 


The  theoretical  development  of  the  reflection/ 


transmission  problem  for  P-SV  waves  is  very  similar  to  the 
development  for  SH  waves.  All  the  equations  and  discussions 


' 
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presented  for  the  SH  case  have  analogies  in  the  P-SV  case. 

It  is  more  complicated  only  in  the  sense  that  there  are  more 
waves  to  handle. 

Since  the  treatment  of  the  P-SV  case  is  analogous  to 
the  SH  case,  we  will  not  expound  it  in  detail.  Instead, 
we  will  confine  ourselves  to  the  derivation  of  the  systems 
of  equations  describing  the  boundary  conditions,  whose 
solutions  are  the  reflection  and  transmission  coefficients 
for  P-SV  plane  waves  impinging  upon  a  plane  boundary 
separating  two  linear  viscoelastic  media.  Free  surface 
reflection  coefficients  will  also  be  treated.  The  notation 
will  be  the  same  as  in  the  previous  section  except,  in 
this  case,  we  add  the  subscripts  P,  (for  a  and  S,  ip  or  8  to 
denote  quantities  associated  with  P  or  S  waves. 

From  equation  (3.42)  we  see  that  the  displacement 
vector  for  the  incident  and  reflected  P  waves  can  be  written 
as 


i  (cot-k 


u 


Pj 


D  .e 
PD 


Pj 


r) 


■pj 


j  =  l,2 


(4.30) 


where,  as  in  (4.1), 


and 


kx  +  (-l)-^dz 

x  a 


d 

a 


+  p.  v. 


(4.31) 

(4.32) 


The  vector  tp^  is  a  complex  unit  vector  representing  the 
particle  motion,  and  is  given  by 


' 
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-*P3 


5?j  kpj  k 


(4.33) 


with 


fcPj 


t„  .  =  1 

p: 


(4.34) 


As  before ,  j— 1  and  2  are  associated  with  the  incident  and 
reflected  waves,  respectively.  The  relative  displacement 


D 


Pj 


Pj 

IS 

Pj 

as 

-i  |  k. 

Pj lBpj 


-ik_B_  . 
P  Pj 


(4.35) 


For  the  transmitted  P  wave,  we  have,  in  the  same  way. 


ul  = 


i  (oot-k  '  »r)  ^ 

D'  e  F  t' 

P  P 


(4.36) 


where 


% 


d' 

a 


H 


k  x 
x 


d 1  z 

a 


+  p . v .  / kl2  -  k2 
—  P  x 

it . 

b  =  _z 

p  ^ 


A 


(4.37) 


%  =  1 


Dp  =  -i  ^ 


The  x-component  of  the  complex  wave  vector  is  the  same  for 
all  waves,  since,  as  before,  we  will  require  that  the 
relative  amplitudes  are  space  independent,  i.e. 


1 


o  z 


kPxl(or 


kSxl>  = 


"Px2 


=  k 


Sx2 


=  k' 
Px 


=  kSx 


=  k 


X 


(4.38) 


As  in  the  previous  section,  this  leads  to  the  complex  form 
of  Snell's  law: 


p^isin0Pi(or 

=  PlsinO'  = 

cp  P 


pn8inesi) 


pJsin0S 


P<j,2sin0P2  =  P<|i2sin0S2 


(4.39) 

Ansin(0pi"YPi)  (or  A<nsin(0srYsi))  =  Vsin(0p2-'rP2) 

=  A<(-2sin<es2‘YS2)  =  A;sin<ep-Yp)  =  A^sinO^-yp 


The  first  equation  can  be  written  in  the  same  familiar 
form  as  (4.14)  by  using  P  =  w/|vj  . 

From  equations  (3.75)  and  (3.76),  noting  that  n=y, 
we  see  that  the  displacement  vector  for  the  incident  and 
reflected  SV  waves  can  be  written  as 


u 


i  (oj  t-k 


Sj 


=  D 


Sj 


Sj 


■  r) 


'Sj 


j  =  l,2 


(4.40) 


where  kc.  =  k  x  +  (-1)  d  z 

o  J  X  p 


dB  =  ±  P-v-^kS  '  kx 


3 

S  x 

•  /\  z\  /\ 

sj 

=  (-1) 

/\ 

yXkSj  '  kSj 

Sj 

•  t  . 

=  1 

S3 

/*v 

*  . 

=  0 

Sj 

S3 

'sj 

=  (-D 

^ i  k_  B  . 

S  S3 

(4.41) 


. 
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For  the  transmitted  S  wave. 


i(mt-k'-r)  ^ 

4  -  4  e  4 


(4.42) 


with  k'  =  k  x  -  d'  z 
S  x  6 


d'  =  +  p.v./ k'2  -  k2 
p  —  S  x 

^  ^  /\  /\  k  * 

4  -  -^xks  -  4  -  kf 


t'  •  t'  =  1 
S  S 


(4.43) 


t*  •  k'  =  ° 


DS  =  "iks  Bs 


The  complex  unit  vectors  can  be  written  in  terms  of  x  and 

/\ 

z  as 


/\  k  /\ 

x 

pj  kp  x 


(-1)  jd 


a 


tp  =  iTT-  X  -  ^  Z 

P  KP  *P 


d' 

a 


3 


(-1)  ■’k. 


(4.44) 


‘Sj 


x  - 


x 


/\  (li  ft  /\  k  /N 

=  _1  x  +  _2£  z 

4  ^  X  ^ 


In  the  case  of  perfect  elasticity,  defining  the  unit  vectors 

A 

tpj ,  etc.  (which  are  real  for  perfect  elasticity)  in  such 
a  way  that  their  x-components  are  positive,  makes  the 
reflection/transmission  boundary  condition  equations 
independent  of  the  direction  of  the  z-axis  (e.g.,  upwards 
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or  downwards  in  Figure  1),  i.e.  the  equations  will  look 
the  same  for  both  directions.  This  convention  has  been 

A 

extended  here  to  the  anelastic  case.  For  instance,  t' 

O 

A  A  A  A 

has  been  defined  as  -yxk'  rather  than  +yxk'. 

The  welded  contact  boundary  conditions  imply  that 


u 


1'  u3,  u31  - - 33 


a. 


and  q  _  _  are  continuous  across  the  boundary. 


In  the  boundary  condition  equations,  the  factor 

i  (oot-k  x) 


t 

i(o3t-k*r)i  ^  .  , 

e  | z_q  =  e  x  will  appear  m  every  term, 

and  hence  we  cancel  it  out  from  the  start.  Using  (4.44) 

in  the  displacements  and  stresses,  the  boundary  condition 

equations  can  be  written  in  their  x  and  z  components. 

For  the  x-component,  the  condition  u^=u|  at  z=0 

yields 

ET  (dpi+dp2>  +  ]T  (dsi+ds2>  =  Ef  +  ds  (4'45) 

P  c  F  o 


and  for  the  z-component,  u3=u3  at  z=0  yields 


a 


(DP1-DP2>  +  E7  (DSrDS2>  =  “EJ  °P  +  Ilf  DS 


(4.46) 


The  horizontal  stress  is  given  by 


°31  =  2Me31  =  M(U3,1  +  Ul,3' 


For  u  in  the  general  form 


+  ^  i  (wt-k*  r )  , 

u  =  D  e  t 


we  then  have 


' 
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->  -> 


a 


31 


-iMD[k  t  +  k  t  ]el(a3t  k'r) 
x  z  z  x 


Hence  the  condition  a 


=  a' 


31  31 


at  z=0  leads  to 


-2Mkx(  )(Dpi-Dp2)  -  M(  BX)(Dgl-Ds2) 

Jr  O 

d'  d  o  2 — k  2 

=  -  2M'kx(^)D^  - 


(4.47) 


The  vertical  stress  is  given  by 


°33  =  A0  +  2Me33  =  A<u1,1+u3,3>  +  2Mu3,3 


For  u  in  the  general  form 


i  (cot-k  •  r )  2 
u  =  De  t 


we  then  have , 


o 


33 


-i { AD  [k  t  +k  t  ]  +  2MDk  t  le1 (wt  k 


x  x 


z  z 


z  z 


r) 


Using  the  relation 

A  =  pa2-2M  =  (M/32)a2  -  2M 
=  M (a/3) 2  -  2M 
=  M  (ks/kp) 2  -  2M 

to  eliminate  A,  the  condition  at  z=0  gives, 

after  some  calculation, 


. 


' 
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2  2 
ci  ^ 

(DP1+DP2>  -  2lVkJ>  <DS1+DS2> 
,  2  2 

=  M,(-n^>DP  -  2M'ds<Jf>Ds 


(4.48) 


Equations  ( 4 . 45) -  ( 4 . 48 )  are  the  equations  describing  the 
P-SV  boundary  conditions.  They  can  be  split  into  two 
systems  of  4  linear  equations  in  4  unknowns,  for  the  case 
of  incident  P  waves,  and  the  case  of  incident  S  waves. 

For  incident  P  waves,  we  must  set  Ds^=0.  The  resulting 
4x4  system  can  then  be  solved  for  the  reflection  and 
transmission  coefficients,  i.e.  (Dp2/D  p)  ,  (Ds2//DPl^  ' 
(Dp/Dp^)  and  (Dg/Dpp) .  Similarly,  for  incident  S  waves, 
we  must  set  Dpp=0.  The  resulting  4x4  system  can  then  be 

solved  for  (Ds2/Dsi^  '  ^Dp2//DSl^  '  ^DP//DS1^  and  ^DS^DS1^  ’ 

The  reflection  coefficients  for  a  free  surface  are 

obtained  from  the  condition  of  vanishing  stress  at  the 
surface.  For  a  P  wave  incident  on  a  free  surface  (D  ^  =  0 )  , 
equations  (4.47)  and  (4.48)  give 


2  2 

2k  d  d^-k 

x  a  fp  -P  )  -  (— — )  n  =  n 

~Y~~  (  PI  p2  {  k^  ’  US2 


d2-k2  2d  k 

<4^>  (DP1+DP2)  -  °S2  ■  0 


(4.49) 


The  solution  of  this  system  gives  the  reflection  coefficients 
for  a  P  wave  incident  on  a  free  surface,  and  they  are 


i , 
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D 


P2 


D 


PI 


D 


S2 


D 


PI 


4kxdadg  -  (drkx)2 

g(kx) 


4kxWkP >(dB-kx> 

g(kx) 


(4.50) 


where  g  (kx)  =  4k2dadg  +  (d2-k2) 2 


(4.51) 


For  an  S  wave  incident  on  a  free  surface  (Dp^=0) , 
equations  (4.47)  and  (4.48)  give 


2k  d  d?-k2 

h4^)DP2  -  <4^  (DS1-DS2)  =  0 


dg-k;2 

°P2 


2dRk 
p  x 


(DS1+DS2)  =  0 


(4.52) 


Solving  this  system  for  the  reflection  coefficients  gives 


D 


P2 


D 


SI 


4kxVkP/kS>  (dB-kx> 


(4.53) 


D 


S2 


D 


SI 


2  2  2  2 

-4kxdad6  +  <Vkx> 


where  g(k  )  is  given  by  (4.51). 

As  mentioned  in  the  previous  section,  one  must  choose 
the  signs  of  da ,  d^,  d^  and  d^  properly  if  one  is  computing 
reflection/transmission  coefficients  in  super-critical 


regions 
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CHAPTER  5 


ASYMPTOTIC  RAY  THEORY  FOR  LINEAR  VISCOELASTIC  MEDIA 


In  this  chapter,  we  apply  asymptotic  ray  theory  to 

the  problem  of  body-wave  propagation  in  a  layered  homogeneous 

isotropic  linear  viscoelastic  medium.  Our  main  goal  is  to 

calculate  an  expression  for  the  geometrical  spreading  of  a 

given  viscoelastic  wave  generated  by  a  point  source  on  the 

surface.  We  develop  only  the  relevant  aspects  of  the  theory 

which  are  necessary  for  our  purposes.  The  development  will 

parallel  standard  applications  of  asymptotic  ray  theory  to 

v 

elastic  wave  propagation  (see,  for  example,  Cerveny  and 
Ravindra  (1971)  and  Hron  and  Kanasewich  (1971)  )  . 

The  equation  of  motion  for  linear  viscoelastic  waves 
in  a  homogeneous  isotropic  medium  is  given  by  equation 
(2.6),  which  is  repeated  here  for  convenience,  i.e. 

[X(t)+y  (t)]*d(V(V-u) )+y  (t)*d(V  u)  =  pu  (5.1) 


A  general  solution  of  this  equation  may  be  written  as 


->■  -> 


u  (r  ,  t)  =2 


7T 


\  c  /  \  iw  [t-x  (r ,  w)  ]  ,  ,  c 

W  (r ,  w)  S  (  go  )  e  do)  (5.2) 


where  S(oo)  is  the  spectrum  of  the  source  pulse,  and  t  is 
the  complex  phase  function  whose  real  and  imaginary  parts 
correspond  to  the  arrival  time  and  the  attenuation  of  a 
ray,  respectively.  W  includes  quantities  pertaining  to 
reflection  and  transmission  at  boundaries,  geometrical 
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spreading,  etc.  To  solve  equation  (5.1)  in  the  sense  of 
asymptotic  ray  theory,  we  expand  W  as  a  power  series  m 


(1/ioo)  with  coefficients  Wn(r,oj),  which,  in  general. 


depend  on  frequency  for  viscoelastic  waves.  We  then  have 

(5.3) 


t 


u ( r , t )  =  —Re 

TT 


<0  i 


£  W  (r , oo )  s  (?,oo)  dco 

On  n 


where 


,  r  ,  S  (oo)  loo? 

s  ( ?  ,  oo)  =  - - —  e 

n  n 

(ioo) 


?  =  t  -  t  (r ,oo) 


(5.4) 


Substituting  (5.3)  into  (5.1),  and  using  the  facts 
ds 


s '  = 


n 


n  d?  Sn-1 


Vs  =  -s  Vt 
n  n-1 


(5.5) 


and  multiplying  through  by  tt  ,  we  obtain,  after  some 
calculation , 


P  Re 


-> 


oo  i 


Z  W  s  -.doo  =  Re 
n  n-2 
n 


[A+y]*d{ 


oo  i 


Z  [s„  9 (W  • Vt ) Vt 
n-  2  n 
n 


-  s  ,(V(W  •  Vx )  +  ( V  •  W  )  Vx  )  +  s  V(V-W  )]dco} 
n— 1  n  n  n  n 


+  y  *d  { 


oo  i 


Z  ?W  (Vt)2-S  1 (2VT-VW  +W  V2t) 

n-2  n  n-i  n  n 

n 


+  s  V  W  ]  dco } 

n  n 


The  general  term  in  this  equation  has  the  form 


(5.6) 


4 


i 
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K t)*d{ 


s  (£,w)F  (r,(o)dw}  ,  m=  0,1,2 

n-m  n 


(5.7) 


(JL)  1 


and  can  be  simplified  as  follows: 

t  °° 


£*an  = 


,-y 


sn-m^n_T  ,U)S)  Fn^r  dco^d0 


OJi 


l (t-n )  { 


0)1 


sn-m-l  O-T  f  m)  Fn  (r  ,03)  dm}dr) 


-e  (4>)  { 


s  1  (t-i-x  ,03)  F  (r  ,03)  do3}dcf) 

n-m-1  n 


0 


03i 


{  io3 


(where  cj)=t-n) 

l(  cj> )  e  iaJ^d(j)  }sn_m  (^  ,03)  F^  (r  ,03)  do3 


03  1 


0 


(where  (5.4)  was  used) 


Letting 

L(w)  =  i  03 
as  in  (2.8),  we  then  have 
t ( t ) *d{ 


t  ( t )  e  "''a)tdt 


0 


s  (£  ,03)  F  (r,o3)do3} 
n-m  ^  n 


03  1 


s  (^,03)L(03)F  (r,co)do3 
n-m  n 


(5.8) 


031 


Hence,  using  (5.8)  and  (2.8),  equation  (5.6)  becomes 
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Re  £{ 
n  ' 


Sn-2  (Pwn)  daJ^  =  Re  Z  { 


00 

f 


-> 


0)1 


n  s 

0)  ! 


sn_2  (A+M)  (Wn*VT)Vxdo) 


sn_1(A+M)  [V  (Wn«  Vt)  +  (V*W  )  Vt]  do) 


+ 


0)j 

CO 


0)1 


s  (A+M)  V  (V  •  W  )  do)  + 
n  n 


s  MW  (Vx)2do) 
ti  z  n 


0)j 


0)1 


S  , M  ( 2Vx •  VW  +W  V2x)do)  + 
n-1  n  n 


s  MV2W  do)} 
n  n 


0)i 


Gathering  terms,  we  can  write  this  equation  as 


00 

Re  I 
n=0 


00 

r 


0)! 


[Sn-2N(T/Wn}  "  sn_i®(T^n)+snL(Wn)  ]do)  =  0 


(5.9) 


where 


N(x,Wn)  =  (A+M) (Wn«Vx) Vx  +  [M(Vx)2-p]Wn 


M(x,Wn)  =  (A+M) [ (V • W  ) Vx  +  V(W  «Vx)] 


+  M[2(Vx*V)W  +  W  V2x] 
n  n 


L(W  )  =  (A+M) V ( V • W  )  +  MV2W 
n  n  n 


(5.10) 


Expanding  the  sum  in  (5.9)  and  defining  W_^  =  W_^  =  0,  we 
see  that  (5.9)  can  be  written  as 


GO 

Re  I 
n=-2 


s  (£  ,oj)  [N  (x  ,W  9)  -M  (x  ,W  )  +L  (W  )  ]  do)  =  0  (5.11) 

n  n+z  n+l  n 


0)i 


In  order  for  (5.11)  to  be  generally  true,  we  must  then  have 


^(T'^n+2}  "  ^(T'^n+l}  +  ^n}  =  0  '  n=-2,-l,...  (5.12) 


i 
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or, 

N(t,Wq)  =  0 

N(t,Wx)  -  M(t,Wq)  =  0  (5.13) 

N(T,Wn+2)  "  M(t,  Wn+1)  +  L(Wn)  =  0,  n=  0 , 1 , 2 , . . . 


Equations  (5.13)  could  have  been  more  easily  obtained 
by  substituting  a  Fourier-transformed  displacement  solution, 
such  as 


-> 

u 


"  +  +  -iu)T(r,w) 

Z  W  (r,u)  2 - - - 

n=0  (100) 


S  (CD) 


into  the  Fourier-transformed  equation  of  motion  (2.7), 
however,  it  was  felt  that  the  above  approach  was  more 
generally  valid. 

The  first  equation  in  (5.13)  implies 

O  O 

N(t,W0).Vt  =  (A+M)  (WQ • Vi )  ( Vt )  +  [M  (Vt)  -p]  (Wq-Vt) 

=  [  (A+2M)  (Vt) 2-p]  (W0-Vt)  =  0  (5.14) 

and  N(t,Wq)xVt  =  [M(Vt)  -p](WqxVt)  =  0  (5.15) 


In  general,  Wq*Vt  and  W^xVt  are  not  zero  simultaneously. 

2  2 

Similarly,  (A+2M) (Vt)  -p  and  M(Vx)  -p  are  not  simultaneously 
zero,  hence  (5.14)  and  (5.15)  have  two  solutions,  i.e. 


(  Vt  )  = 


P  _  1 


A+2M 


a 


WqxVt  =  0 


(5.16) 


(Vt) 


2  _  p  _  1 


M 


f3‘ 


Wq»Vt  = 


r 


0 
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These  are  the  eikonal  equations.  In  terms  of  previous 
notation,  x  is  given  by 


k»r  _  P • r  _  .  A»r 

00  00  CO 


(5.17) 


for  either  P  or  S  waves,  hence 


->  ->  -> 

„  k  P  .  A 
Vx  =  —  = - i  - 

co  oo  co 


(5.18) 


=>  (Vx)  = 


k  »k 


co 


k2 


and  so  equations  (5.16)  are  the  same  as  equations  (2.11). 


The  amplitude  coefficients  W^  can  be  extracted  from 
equations  (5.13).  Since  we  will  be  taking  a  plane  wave 
approach,  it  will  not  be  necessary  to  determine  for 
all  n,  rather  only  for  n=0,  i.e.  we  consider  only  the 


zeroth  order  of  the  theory.  We  now  calculate  Wq  for  P  waves 


5.1  P  Waves 


For  plane  P  waves,  we  can  write 


W0  =  w0tp 


(5.19) 


where  t  is  the  complex  unit  vector  for  P  waves,  discussed 
in  section  4.2,  and  is  given  by 


tp  =  kp  =  aVx 


To  determine  WQ,  we  first  compute  M(t,Wq)*Vt 


(5.20) 


X 
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M(t,W0).Vt  =  (A+M) [ (V*WQ) (Vt)  +  V ( WQ • Vt ) •  Vt  ] 

2  ^  ^ 

+  M [V  t(Wq»Vt)  +  2 { ( Vt • V ) Wq } • Vt ] 

We  have 

(Vt)2  =  1/a2 

V*WQ  =  VMW^aVT)  =  aVT*VW^  +  aW^V2T 

Wq*Vt  =  W^a(VT)2  =  W^/a 

( Vt • V) WQ  =  ( Vt • V ) WQtp  =  (VT*VWQ)tp 

(since  the  medium  is  homogeneous) 
=>  Vt  •  {  (  Vt  •  V )  Wq  }  =  (1/oi)Vt*VWq 


Hence 

M(t,WJ -Vt  =  (A+M)[— Vt*  VWn+— Wn V^t+— Vt • VWn ] 
0  a  0  a  0  a  0 

W 

+  M  [— V2t  +— Vt  •  VWj 

a  a  0 

=  (//A  [WgV2!  +  2Vx  •  VWQ ] 

=  pc([WqV2t  +  2Vx*VW„] 


^  ^  2 
Hence  [M (t ,Wq) • Vt] =  paV* [ (W^)  Vt] 

Now,  from  (5.10),  we  have 

N(T,Wn)*VT  =  (A+M)  (Wn*VT)  (Vt) 2  +  [M  ( Vt ) 2-p ]  (Wr • Vt ) 


r  ,  A  +  M  \  ,  M  1  //r  V7  \ 

=  [<  — 2>  +  -2  -Pi  <Wn-VT) 

a  a 


r  A+2M 
1 — 

a 


p]  (W  *Vt) 
n 


0 


(5.21) 

(5.22) 


(5.23) 


1 
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Hence,  the  second  equation  of  (5.13)  then  shows  that 
M(t,Wq)*Vt  =  0.  Therefore,  equation  (5.22)  becomes 

V* [ (WQ) 2Vt]  =  0  (5.24) 

Let  us  now  use  the  concept  of  a  ray  tube  about  a 
central  ray  connecting  two  wavefronts  of  the  central  ray 
at  different  times  tQ  and  t  with  t  <t  (see,  for  example, 
Cerveny  and  Ravindra  (1971) ) .  Let  s  be  the  length  along 
the  ray  path  in  the  direction  of  propagation  P.  The  ray 
tube  is  bounded  by  rays  adjacent  to  the  central  ray  and 
infinitesimally  distant  from  it,  and  by  the  infinitesimal 
wavefront  surfaces  at  s=sq  and  s,  i.e.,  da  (sq)  and  da  (s) . 
Integrating  over  the  volume  V  of  the  ray  tube,  and  using 
Gauss'  divergence  theorem,  we  have 

2 

(WQ)  Vt  *da  =  0  (5.25) 

y 

a 

where  da  is  the  element  of  area  on  the  ray  tube.  Writing 
Vt  as 

■j  ^  i  A  A 

Vx  =  -  (P-iA)  =  -  (PP-iAA)  (5.26) 

03  03 

/\  /S 

where  P  and  A  are  real  unit  vectors,  we  obtain 

(Wn)2Vx*da  =  —  { [ (Wn)  da]  -  [ (Wn)  da]  } 

0  03  U  s  u  =>0 

.  i^£iT{[(w0)2da]s  -  [(w0)2da]So3 


0 


(5.27) 


h 
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Hence 

[ (WQ) 2da ] s  -  [(WQ)2da]s  =  0 

o 


which  leads  to 


W0(s)  = 


da  (s  ) 
o 

da  ( s ) 


*\ 


Vso} 


(5.28) 


(5.29) 


This  equation  gives  at  any  point  s  on  the  ray,  assuming 

we  know  it  at  s  .  It  has  the  same  form  as  the  analogous 

o  ^ 

.  v  ^ 

elastic  result  (see  Cerveny  and  Ravindra  (1971)),  which  is 
what  one  would  expect  since  (5.29)  is  a  purely  geometrical 
result  which  describes  the  geometrical  spreading  of 
waves.  Intuitively,  one  would  also  expect  an  attenuation 
term  to  be  present,  however,  this  is  accounted  for  by  the 
imaginary  part  of  x . 


5.2  S  Waves 

For  S  waves,  we  can  write 


2  -*  3  -> 

W  t  +  W  t 
w0  2  0  3 


(5.30) 


where  the  complex  unit  vectors  t ^  an<3  t^  are  given  by 

t2  =  yxks  ,  t3  =  y  (5.31) 

as  discussed  in  section  4.2.  The  vectors  t  ,  t ^  and  t^ 
are  orthogonal.  To  determine  WQ,  we  compute  M(T,WQ)«t^, 
i=2 , 3 .  Since  ti«Vx=0,  and  Wq*Vt=0,  we  have 

M(x,W0)-ti  =  M[WqV2t  +  2 { ( Vx • V) WQ } • ti ]  ,  i=2 , 3 


■ 


1 
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Now, 

{ (VT-V)W0}*ti  =  (Vt-VWq) t2*ti  +  (Vt-VWq) t  *t 

=  Vt • VWq  ,  i=2 , 3 

Hence 

M(T,W0)*ti  =  M[WqV2t  +  2V  t • VWq ]  (5.32) 

and  [M(t,Wq) -tjwj  =  MV-[(Wq)2Vt]  (5.33) 

2  2 

Now,  since  (Vt)  =  1/3  =  p/M,  we  have  for  S  waves, 

N(x,Wn)  =  (A+M) (W  • Vt ) Vt  (5.34) 

Hence  N(T,51)*ti  =  0  ,  i=2,3  (5.35) 


Therefore,  the  second  equation  of  (5.13)  shows  that 
M (t  , Wq)  •  t_^  =  0.  Hence,  equation  (5.33)  gives 


V- [ (Wq) 2Vl]  =  0 


(5.36) 


As  in  the  previous  section,  this  leads  to 

,  i=2,3 


W0(s)  = 


da  (s  ) 

daTif  W0(so> 


(5.37) 


as  in  equation  (5.29)  .  It  again  has  the  same  form  as  the 

V 

analogous  elastic  result  (see  Cerveny  and  Ravindra  (1971) ) . 


5.3  Geometrical  Spreading 

We  now  calculate  the  geometrical  spreading  factor  for 
a  ray  propagating  through  a  layered  homogeneous  isotropic 
linear  viscoelastic  medium.  Let  us  consider  a  typical  ray 


. 


(t 
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consisting  of  m  segments,  as  in  Figure  8  and  define  0.  to 
be  the  endpoint  of  the  jth  ray  segment,  at  an  interface. 
Using  equations  (5.29)  or  (5.37)  to  compute  the  amplitude 
coefficient  of  the  ray  anywhere  along  its  path,  we  obtain 
at  the  endpoint  K 


W  (K) 


1 

L(K,Kq) 


m-1 

n 

j=i 


R  . 
3 


where 


L(K,Kq) 


_ /do  (K) 

da  (K  ) 

o 


m-1 

n 

j=i 


'da (0  . ) 
da ' (0^ ) 


(5.38) 


(5.39) 


We  assume  that  K  is  at  a  unit  distance  from  the  source. 

o 

R.  is  the  reflection/transmission  coefficient  at  0 . .  We 
3  3 

use  the  prime  to  denote  quantities  associated  with  the 
reflected  or  transmitted  wave  at  0_.  ,  and  an  unprimed 
quantity  is  associated  with  the  incident  wave  at  Cu  . 

L(K,K  )  is  the  geometrical  spreading  factor  of  the  ray. 

To  evaluate  (5.39),  we  first  get  an  expression  for 
da (K) /da (K  ) .  Referring  to  Figure  9,  we  see  that 


da  (K)  =  (cosG  (K)dx.)  (xd(J)Q) 


where  d)  is  the  azimuthal  coordinate  at  K  ,  and  x  is  the 
o  o 

epicentral  distance.  Now  since  x=x  ( z  ,'0  ,  y  )  where 

0  =0 (K  )  and  y  =y(K  )  are  the  propagation  and  attenuation 
o  o  1  o  '  o 

angles  at  the  source,  we  have 


dx  = 


3x 

30 


d0 


+ 


3x 

W, 


dy 


3x 

sir 


dz 


■ 


n 
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dZ 


Oj 


FIGURE  8:  Typical  ray  in  a  homogeneous  layered  medium. 


K  and  K  are  the  source  and  receiver  points , 
o 

respectively.  CK  ,  j  =  1,  ...,  m  is  the  endpoint 

of  the  jth  ray  segment,  at  an  interface. 


- 
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FIGURE  9:  Ray  tube  between  Kq  and  K  for  a  homogeneous 

medium.  da(K)  is  the  cross-sectional  area  at 
K.  The  diagram  is  partially  schematic  -  there 
may  be,  in  actuality,  an  arbitrary  number  of 
reflections  and  transmissions  at  interfaces 
between  K  and  K. 
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From  Figure  9,  dz=0.  Also,  we  will  assume  that  the  source 

rays  all  have  the  same  value  of  yq,  which  implies  dyo=0. 

The  cross-section  of  the  spherical  wave  at  the  source, 

da  =da(K  ),  is  given  by  da  =sin0  d0  dd>  .  Hence,  we  have 
o  o  ^  1  o  ooyo 


da  (K) 

da  (K  ) 

o 


x(3x/30  )cos0(K)d0  deb 
o  o  o 

sin0  d9  3/ 

o  o  o 


x ( 3x/30q) cos0 (K) 
sin0 

o 


(5.40) 

Next,  we  simplify  the  product  in  (5.39).  From  Figure  10, 
we  see  that 


da (0 . )  cos0 (0 . ) 

_ 3 _  =  3 

da'  (0.)  cos0 '  (0 . ) 
3  3 


(5.41) 


Now,  for  homogeneous  media,  0  1  (0 ^  )  =0  (0 _.+^)  •  Using  this  in 
(5.41)  and  then  expanding  the  product  in  (5.39) ,  we  obtain 


m-1  /da (0 . ) 

n  ./  ] 


j-i 


cos0  (0-^) 

cos0 (0  T 

m 


COS0 
_ o 

cos0 (K) 


(5.42) 


Hence,  using  (5.40)  and  (5.42),  equation  (5.39)  becomes 

=  /x  cot9  (5.43) 


L(K,Kq) 


30  COt0o 


The  epicentral  distance  x  is  given  by 


m 

x  =  Z 

j  =  l 


x 


(5.44) 


where  x  .  =  h .  tan  0 . 

3  3  3 


(5.45) 


where  h.  is  the  thickness  of  the 
3 

ray  segment,  0^  =  0  (Ch  ) ,  and  x^  is 
traversed  by  the  jth  ray  segment. 


layer  containing  the  jth 
the  horizontal  distance 


. 


I 


FIGURE  10:  Change  in  the  cross-sectional  area  of  the 


ray  tube  at  an  interface. 


' 
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We  now  evaluate  3x./90  • 

3  o 


tan0  . 
3 


P  . 
3* 


P  . 
jz 


We  have 


where 


P  .  sin0  . 
3  3 


P  sin0 
o  o 


d.  =  +  p.v./k^-k^ 
3  3  x 


(P ,  SV  or  SH  waves) 


(5.46) 


(5.47) 


We  have  suppressed  the  subscripts  P,  S,  a, 6  etc.  so  that 
the  analysis  holds  for  P,  SV  and  SH  waves.  The  subscript 
"o"  denotes  quantities  at  K  .  The  "+"  sign  was  chosen  for 
dj  since  we  are  dealing  only  with  body  waves,  i.e.  ray 
segments  propagating  at  sub-critical  angles.  Hence,  we 
have 

x.  =  h.P  sine  /d._  (5.48) 

3  3  °  o  jR 

which  implies 


'x 


J  =  h.P  [_2R  ° 

9  0  j  o 


d.^cos0  -sin0  (9d._/90  ) 


o 


o  jR  o' 


-] 


(5.49) 


For  (9d.  /3 0  )  we  have,  suppressing  the  symbol  "p.v.", 
jR  o 


9d  k  9k 

.  JR  =  Re  r-  __*] 
90  Ke  L  d.  90  J 

o  j  o 


(5.50) 


Now  k  and  9k  /30  are  given  by 
x  x  o 


k  =  P  sin0  -  iA  sin(0^-y^) 
x  o  o  o  o  o 


3k 

=  P  cos0  -  iA  cos(0  -y  ) 
90  o  o  o  o  o 


(5.51) 
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Substituting  these  into  (5.50),  and  then  using  the 
resulting  equation  for  3d_.R/30Q  in  (5.49),  we  obtain 


3x  .  h  .P 

— — {cos0  + 
30  d  o 

o  jR 


2  2 

h[~P  sin20  -A  sin2(0  -y  )  ]  sin0 
o  oo  o  o  o 

l^i 

1  D 


P  A  sm  20  -y  )sm0  d._ 
o  o  o  'o  o  jIt. 

I  25  J 

d.D  dz 

0R1  : 1 


(5.52) 


Hence,  we  can  write 


3x 


m  3x  . 

=  E  =  (P  cos0  )  Z 

.  .30  o  o 

3  =  1  o 


(5.53) 


where 

m  h .  ~  9  m  h  . 

Z  =  E  -r-3—  +  ^[PZsin20  -A~sin2(0  -y  )  ]  tan0  E  - 3 


j=l  djR 


o  o 


o  o 


°  j=l  d.  |dr 
J  jR1  : 


m  h  .d 

-  P  A  sin  (20  -y  )  tan0  E  — 2 — 3 — 
°  o  o  o  °  i=l  d2  Id2 

J  jR'  j 


m 


Using  x  =  E  x.  =  P  sin0 
j  =  l  D 


m  h . 

E  J- 


j=l  djR 


(5.54) 


(5.55) 


implies 


m  h . 

xcot0  =  P  cos0  E  — 3- 
o  o  o  d^ 


(5.56) 


j=l  jR 

Hence,  the  final  expression  for  the  geometrical  spreading 


factor  L  is 


L(K,Kq) 


P  COS0 
o  o 


m 
(  E 


h  . 


•  i  d  .  I, 

j  =  l  3R 


3-)  Z 


(5.57) 


where  Z  is  given  by  (5.54) .  We  can  show  that  this  reduces 
to  the  correct  expression  for  L  in  the  perfectly  elastic 


■ 


105 


case,  as  follows.  For  perfect  elasticity,  we  have 


d  .  =  d .  =  —  cos0 . 

3R  5  vj  : 


A  =  0 
o 


Hence  (5.54)  becomes 


2  2 

1  m  v . h .  ,  sin  0  m  v  h . 

7  =  —  y  J  5  _L  ±  _ o  y  j  j _ 

oo  cos0  .w  2  3 

j-1  3  v  3=1  cos  0 . 


m  v  .  h  . 


—  £  [—5 — - — ^ —  (v  cos  0.  +  v.sm  0  )] 

CO.,  2  3,  O  1  1  O 

3  =  1  v  cos  0 .  J  J 

o  j 


where  Pq=co/vo  was  used.  Now,  Snell's  law  in  elastic  media 

2  2 

implies  v.sm0  =v  sin0  .  ,  hence,  using  cos  0 .+sin  0.  =  1, 
j  o  o  j  y  3  j 

Z  becomes 

,  m  v  .  h  . 

Z„,  =  i  E  — 2  J 
et  co  .  ,  3n 

3=1  cos  0 . 

J 


and  finally,  (5.57)  becomes 

cos0  / 


'<Ll 


v 


o 


m  v  .  h  .  m  v  .  h  . 

(  £  ^2-)  (  z  -U-) 

3=1  3  3=1  cos  0 . 

3 


(5.58) 


This  agrees  with  the  expression  for  L  in  perfectly  elastic 
media  given  by  £erveny  and  Ravindra  (1971,  equation  2.143) 
The  expression  (5.57)  for  L  will  be  used  in  the 
next  chapter  for  the  computation  of  body-wave  synthetic 
seismograms  for  a  point  source  at  the  surface. 


i  ' 


> 


CHAPTER  6 


SYNTHETIC  SEISMOGRAMS 

In  this  chapter,  we  use  the  theory  developed  in  the 
previous  chapters  to  compute  synthetic  seismograms  for  body 
waves  in  a  layered  anelastic  medium.  Both  the  teleseismic 
case,  and  the  case  of  a  point  source  at  the  surface  are 
treated.  Hron  et  al .  (1974)  developed  a  computer  program 

to  calculate  synthetic  seismograms  for  teleseismic  body 
waves  in  a  plane  layered  homogeneous  isotropic  elastic 
medium  via  a  ray  theory  approach.  A  similar  program  to 
compute  synthetic  seismograms  for  waves  generated  by  a 
point  source  at  the  surface  was  developed  by  Hron  and 
Kanasewich  (1971) .  These  programs  automatically  generate 
all  rays  which  arrive  at  the  receiver,  including  multiples, 
up  to  a  given  maximum  number  of  segments  per  ray,  and  can 
automatically  identify  the  rays  that  make  up  a  given  event 
on  the  synthetic  seismogram.  The  program  for  the  teleseis¬ 
mic  case  can  also  handle  dipping  interfaces.  The  program 
for  the  surface  point  source  case  groups  all  generated  rays 
into  kinematic  and  dynamic  analogs,  which  significantly 
improves  the  computing  efficiency.  For  the  purposes  of 
this  thesis,  these  programs  have  been  modified  to  compute 
synthetic  seismograms  for  SH  body  waves  propagating  through 
a  layered  homogeneous  isotropic  linear  viscoelastic  medium. 
In  order  to  be  able  to  compute  synthetic  seismograms  for 
a  layered  medium  consisting  of  elastic  layers  embedded 
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among  anelastic  layers,  all  six  computational  cases  of  re¬ 
flection/transmission  at  a  boundary  discussed  in  Chapter 
4  have  been  incorporated  into  the  programs. 

The  medium  we  use  is  a  model  of  the  Berkeley  crust 
due  to  Silva  (1976),  and  is  shown  in  Table  1.  The 
parameters  v  and  vuo  are  the  homogeneous  phase  speeds 
for  P  and  S  waves  and  the  parameter  h  is  the  thickness  of  a 
layer.  The  model  consists  of  three  anelastic  layers  over 
an  elastic  half-space.  The  interfaces  are  plane  and 
horizontal . 


6 . 1  Teleseismic  Case 

For  teleseismic  waves,  each  ray  originates  in  the 
elastic  half-space,  which  means  that  A  will  be  zero  for  the 

X 

incident  ray  segment,  and  since  A  must  be  continuous  across 

X 

a  boundary,  A  will  be  zero  for  every  ray  segment,  i.e.  A 

X 

will  be  in  the  +  z  direction  for  every  ray  segment. 

For  SH  waves,  the  displacement  of  a  given  ray  at  the 
receiver  is  in  the  y  direction,  and  the  amplitude  of  motion 
of  the  ray  is  given  by  an  inverse  Fourier  transform  as 


W  = 


-  Re 

TF 


,  ,  i c()  ( oj )  _  ,  «  iw  ( t-T  (go)  ) 

Y  (go)  e  S  (go)  e  doo 


(6.1) 


=  -  Re 

TT 


.  GOT  ico(t-T) 

Ye  9S(G0)e  e  dGo 


o 

oo 

We  may  use  the  above  form  rather  than  (1/2tt)  /  because  W 

—  OO 

is  the  response  and  is  a  real  number  (see  Appendix  2) . 
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TABLE  1 

CRUSTAL  MODEL 


Layer 

Number 

VHP 

(km/s) 

VHS 
(km/s ) 

1 

4.2 

2.4 

2 

6 . 1 

3.5 

3 

7.3 

4.2 

4 

7 . 8 

4.5 

P 

(g/cc) 

Qp 

Qs 

h 

(km) 

2 . 1 

67 

30 

1.4 

2.6 

100 

45 

8.2 

3.0 

180 

80 

12.9 

3.3 

CO 

CO 

00 

. 

' 


is  the  product  of  the  anelastic 


The  quantity  Y  (go)  e^  ^ 
reflection  and  transmission  coefficients  for  the  ray,  Y  and 
(j)  being  the  relative  magnitude  and  phase.  S  ( oo )  is  the 
spectrum  of  the  source  pulse.  t  =  x  +ix  is  the  complex 
phase  function,  where  x  is  its  real  part  and  is  the  arrival 
time  of  the  ray,  and  x^.  is  its  imaginary  part  and  correspond 
to  the  attenuation  of  the  ray.  For  a  given  ray  segment, 
x  is  given  by 


->■  ->  -*• 

k  *  r  _  P • r  .  A*  r 

CO  ~  03  1  _ CO 


(6.2) 


hence , 


x  =  P  •  r/ co ,  t  =  -A*  r/co  .  (6.3) 

R  X 

We  suppress  the  subscript  "S",  since  we  are  dealing  only 
with  SH  waves. 

Since  we  are  considering  teleseismic  waves,  no 
geometrical  spreading  factor  appears  in  (6.1). 

In  our  computations,  we  assume  that  Q  is  independent 
of  frequency.  This  means  that  we  are  neglecting  velocity 
dispersion  in  our  calculations,  i.e.  we  are  using  one 
constant  velocity  v  for  all  frequencies.  It  is  well 
known  that  this  assumption  introduces  an  acausality  in  the 
response  (Futterman  (1962)).  However,  this  will  not  be 
an  important  effect  in  our  case,  due  to  the  nature  of  our 
source  pulse,  whose  spectrum  exhibits  a  large,  narrow  peak 
at  a  dominant  frequency  f  .  Thus,  most  of  the  contribution 
to  the  integral  in  (6.1)  will  come  from  a  narrow  frequency 
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interval  centered  about  w  =  co  =  2irf  ,  and  the  dispersion 

o  o 

in  this  interval  will  be  negligible.  This  also  means  that 
choosing  our  constant  velocity  v  to  be  the  velocity  at  the 
dominant  frequency  of  our  source  pulse,  i.e.  v  =  v  (u)  ) 
(which  can  be  obtained  from  dispersion  data) ,  and  similarly 
choosing  Q  =  Q(w  ),  should  considerably  reduce  the  errors 
in  arrival  time  that  would  be  introduced  by  neglecting 
dispersion.  Finally,  if  we  restrict  ourselves  to  values  of 
Q  which  are  not  too  low,  and  ray  paths  which  are  not  too 
long,  we  should,  by  the  above  frequency-independent  approach, 
be  able  to  obtain  a  good  approximation  to  the  actual 
dispersive  response.  These  arguments  would  not  be  as 
valid  in  the  case  of  a  6-function  source  pulse,  whose 
spectrum  is  a  constant  for  all  frequencies. 

An  added  computational  advantage  in  using  the  above 
frequency-independent  approach,  is  that  by  choosing  an 
appropriate  source  pulse,  we  will  be  able  to  calculate  the 
response  in  (6.1)  analytically,  i.e.  in  closed  form,  as  we 
shall  see  below. 

The  arrival  time  t  is  obtained  by  summing  up  the 
travel  times  of  each  segment  of  the  ray.  The  speed  and 
length  of  each  ray  segment  are  calculated  from  the 
viscoelastic  theory  presented  in  the  previous  chapters. 

The  attenuation  term  t  is  obtained  in  a  similar 
fashion.  For  a  given  ray  segment,  from  (6.3)  we  have 

VTj  =  -A/o) 


(6.4) 


. 


' 


» 


V 


Ill 


Letting  s  be  the  length  along  the  ray  path  in  the  direction 
of  propagation  P,  and  using  (2.16)  and  (6.4)  we  then  obtain 
for  a  given  segment 


__  =  p  •  Vx x  =  |  Vt I  |  cos  (P  ,  Vt x ) 

A  A 

=  -  cos(TT-y)  =  -  -  cos  Y 

Im(k2) 

2coP 


(6.5) 


Hence,  for  the  entire  ray,  t  is  given  by  a  sum  over  the 
ray  segments: 


t ^  (receiver) 


£ 

j 


rIm(k2) 
L  2 cop 


]  .  ds 
] 


(6.6) 


where  the  limits  of  the  integral  are  the  beginning  and  end 
points  of  the  j-th  ray  segment.  Using  (2.15),  and  the  fact 
that  the  medium  is  homogeneous,  (6.6)  becomes 

x  (rec.)  =  E  [i  (k2/u2)  T  (v.l  )  1  (6.7) 

j  ^  J  J  J 

where  v^  =  co/P^  is  the  phase  speed  of  the  j-th  ray  segment, 
and  is  its  length.  Employing  (2.21),  equation  (6.7) 
finally  becomes 

Q^v.l 

x  (rec.)  =  -  £  [ - - - ^—2 - 1  (6.8) 

j  2  / 

VHj(1+/l+Qj  > 

If  Qj  is  frequency-independent,  then  so  are  v^  and  v^ , 
hence  x  is  frequency-independent.  Also,  Y  (co)  and  4>  (to) 


. 


(6.1)  becomes 


are  frequency-independent.  Letting  E,  =  t-x  , 

R 


W  =  -  Re{e1(J) 

7T 


COX. 


0  ,  s  I  io)£  , 

S  (co)  e  e  dco  j 


o 


=  Y  Re{el4)f  (^xI)  } 

=  Y  [fR  (£  ,  x  I )  cosij)  -  fx  (CfTj)  sin({)] 


(6.9) 


where 

f (?,ti) 


COX 

S  (co)  e 


I  i  co  E, 

e 


dco 


(6.10) 


o 

The  subscripts  R  and  I  refer  to  the  real  and  imaginary  parts. 

It  would  be  a  definite  computational  advantage  to  be 
able  to  evaluate  (6.10)  analytically,  i.e.  in  closed  form, 
since  then  the  response  could  also  be  evaluated  in  closed 
form.  With  this  in  mind,  we  choose  our  source  pulse  to  be 

sin(w  t) 

s(t)  =  - 2 -  (6.11) 

1+  (o)ot/n) 

The  spectrum  for  this  pulse  is  given  by 


S  (co) 


4^-{exp  [-—  I  co-co  |]  -  exp[-~(co+oo  )]} 

2oo  co  1  o'  co  o 


(6.12) 


The  pulse  and  spectrum  are  shown  in  Figure  11.  As  required, 
the  spectrum  has  a  peak  at  co=coq.  Since  x^.  and  E,  are 
independent  of  co ,  we  can  see  that  the  integrand  in  (6.10) 
involves  a  product  of  exponential  and  trigonometric  func¬ 
tions  whose  arguments  are  linear  in  co ,  which  is  easy  to 
integrate.  The  result  is 


. 
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FIGURE  11: 


Source  pulse,  and  absolute  value  of  spectrum 


divided  by  C,  where  C  =  7Tri/2coo. 
diagram,  r)  =  3  and  f  = 
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t^{e  °  ^  r  r  ,  r _  ^-_-rl 


fR(5'TI)=- 


Q  ntQ+£  _TI)  sinw0^+2TI^cosooo^] -2TI^e  } 

[(W2H2]  K^-X^^+C2] 


(6.13) 

03  T 

f  ,P  T  ,_.  o  °  I[(t2+C2-T2)cos(J)oC-2TI5sincoo?]-(t2+?2-T2)e'n} 

r  T  '  ^  '  L  T  '  "  ‘  T T - - 7— - 

[(Vti}  +n  [(t0-Tz)  +ri 

where 


t  =  n/co 

o  o 


(6.14) 


For  x  =0,  i.e.  no  attenuation,  these  reduce  to 


fR(C,°)  = 


sinco  £ 
o 


1+  (cjoC/ti) 


fI(C,0)  =  - 


(cosoo  £-e  r|) 

1+  (<i)o£/n)  2 


(6.15) 


As  is  expected,  f  (£,0)  is  the  same  as  (6.11),  where  £ 

K. 

corresponds  to  t. 

The  synthetic  seismograms  for  the  crustal  model  in 

Table  1  are  shown  in  Figure  12.  Figure  12a  shows  the 

results  for  the  case  of  perfect  elasticity  (Q  ^  is  set 

equal  to  zero  in  each  layer) ,  and  Figure  12b  shows  the 

results  for  the  anelastic  case.  The  source  parameters 

used  are  n=3  and  f  =10  Hz.  The  interfaces  in  the  crustal 

o 

model  are  horizontal.  From  the  top  down,  the  traces  in 
Figures  12a  and  12b  correspond  to  teleseismic  waves 
impinging  upon  the  deepest  interface  at  incident  angles  of 
0°,  15°,  30°,  45°,  60°  and  75°.  In  each  figure,  the 
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FIGURE  12:  Synthetic  seismograms  for  teleseismic  SH  body 
waves 

(a)  for  the  perfectly  elastic  crustal  model 
(Q  ■*"  =  0  in  each  layer,  and 

(b)  for  the  anelastic  crustal  model. 

In  each  figure,  from  the  top  down,  the  angles 
of  incidence  for  the  incoming  ray  in  the  half¬ 
space  are  0°,  15°,  30°,  45°,  60°  and  75°.  The 
scales  are  shown  in  the  upper  right-hand 
corners.  The  letters  a,  b,  c,  etc.,  identify 
the  phases  contributing  the  most  to  the  events 
(see  Table  2) .  The  source  parameters  are 

T)  =  3  and  f  =  10  Hz. 

o 
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scale  for  each  trace  is  shown  in  the  upper  right-hand 
corner.  The  letters  a,  b,  c,  etc. ,  identify  the  rays 
contributing  the  most  to  the  events  on  the  trace,  and  are 
listed  in  Table  2.  The  numerical  code  is  such  that  4321 
stands  for  S^S^S^S^,  etc. 

The  traces  in  Figure  12  show  the  effects  of  attenua¬ 
tion:  the  amplitudes  are  reduced,  the  oscillations  are 
smoothed  out,  and  the  waveform  spreads.  Many  of  the  later 
arrivals  which  appear  in  the  elastic  traces  do  not  show  up 
in  the  anelastic  traces,  and  also,  in  a  given  trace,  the 
ratio  of  the  amplitude  of  a  given  later  arrival  to  the 
amplitude  of  the  first  arrival  is  significantly  smaller  in 
the  anelastic  case  than  in  the  elastic  case.  This  is  to  be 
expected  because  the  amplitude  reduction  is  not  only  due  to 
reflection/transmission  losses,  but  also  to  attenuation, 
which  increases  monotonically  with  the  ray  path  length. 


6.2  Surface  Point  Source  Case 


In  this  case,  the  SH  waves  propagating  through  the 
medium  are  generated  by  a  point  source  situated  at  the 
surface,  and  hence,  the  expression  for  the  displacement  of 
a  given  SH  ray  at  the  receiver,  i.e.,  equation  (6.1),  must 
contain  an  extra  factor  L,  representing  the  geometrical 
spreading  of  a  ray  generated  by  a  point  source,  i.e., 

no 

M 


W  =  -  Re  / 


Ye-  ..  ,  WTI  iw(t-TR> 

— -  S  (oo)  e  e  dm 


(6.16) 
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TABLE  2 

RAY  CODES  FOR  THE  TELESEISMIC 
SYNTHETIC  SEISMOGRAMS 


Letter  Numerical 

Code  Code 


a 

b 

c 

d 

e 

f 

g 

h 


3 

k 

1 

m 

n 

o 

P 

q 

r 


4  3  2  1 


4  3  2  1  1  1 


43211111 


4321111111 


432111111111 


4  3  2  2  2  1 


43211221 


4  3  3  3  2  1 


4321111221 


43332111 


432111111221 


4322211221 


43223321 


432112211221 


4321123321 


432111123321 


432111122111 


432112332111 


. 
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where  L  is  given  by  (5.57)  for  a  linear  viscoelastic  medium. 
The  expressions  for  t  ,  t  ,  S  (co)  ,  etc.  are  the  same  as  in 

K  _L 

the  previous  section,  and  we  again  assume  that  Q  is 
independent  of  frequency  (which  implies  that  L  is  also 
frequency-independent) . 

Since  the  initial  ray  segment  is  located  in  the  first 

layer,  which  is  an  anelastic  layer,  we  are  faced  with  the 

choice  of  the  value  of  y  ,  the  attenuation  angle  of  the 

initial  ray  segment.  Intuitively,  one  would  expect  that 

the  source  would  generate  homogeneous  rays,  for  which  yQ=0, 

however,  in  general,  yQ  is  an  arbitrary  parameter.  We 

assume  only  that  each  ray  generated  by  the  source  has  the 

same  value  of  y  .  The  expression  for  the  geometrical 

spreading  factor  L  in  (5.57),  gives  L  for  any  value  of  y  . 

In  the  examples  below,  we  compute  seismograms  for  three 

different  values  of  y  ,  and  examine  their  differences. 

o 

The  synthetic  seismograms  for  the  crustal  model 

of  Table  1  are  shown  in  Figure  13.  Figure  13a  shows  the 

results  for  the  case  of  perfect  elasticity  (Qg1  is  equal 

to  zero  for  each  layer)  and  Figures  13b,  13c  and  13d  show 

the  results  for  the  anelastic  case  for  yo=-60°,  0°  and 

+60°,  respectively.  The  source  parameters  are  n=3  and 

f  =15  Hz.  Traces  for  eight  epicentral  distances  (0-3.5  km) 
o 

are  shown.  In  each  figure  the  scale  is  shown  in  the  upper 
right-hand  corner.  The  letters  a  to  f  identify  the  rays 
contributing  the  most  to  the  events  on  the  traces,  and  are 
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FIGURE  13:  Synthetic  seismograms  for  SH  body  waves 
generated  by  a  surface  point  source 

(a)  for  the  perfectly  elastic  crustal  model 
(Q  ^  =  0  in  each  layer) ,  and 

(b)  -  (d)  for  the  anelastic  crustal  model. 

The  attenuation  angle  yq  of  the  incident  ray 
segment  is  yq  =  -60°,  0°  and  +  60°  in  (b)  ,  (c) 

and  (d) ,  respectively.  Traces  for  eight 
epicentral  distances  from  0  to  3.5  km  are 
shown.  The  scales  are  shown  in  the  upper 
right-hand  corners.  The  letters  a,  b,  c ,  etc., 
identify  the  phases  which  contribute  the  most 
to  the  events  (see  Table  3) .  The  source 
parameters  are  r\  =  3  and  f  = 
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listed  in  Table  3.  The  numerical  code  is  the  same  as  for 
Table  2,  i.e.,  for  example,  12  2  1  stands  for  S]_S2S2S1* 

In  these  seismograms,  we  do  not  see  as  many 
significant  arrivals  beyond  the  initial  arrival  as  we  did 
in  the  teleseismic  case,  due  to  the  presence  of  geometrical 
spreading,  which  reduces  the  amplitudes  of  the  rays.  As 
for  the  teleseismic  case,  we  can  see  the  effects  of 
attenuation  in  the  reduction  of  amplitudes,  and,  for  the 
rays  with  longer  ray  paths,  the  spreading  of  the  waveform 
and  the  smoothing-out  of  the  oscillations. 

The  arrival  labelled  "a"  in  the  four  figures  is  the 

primary  reflection  off  the  first  interface,  i.e.  the 

ray,  and  the  behaviour  of  its  amplitude  and  phase  for 

increasing  epicentral  distances  illustrates  the  behaviour 

of  the  HlHl  reflection  coefficients  of  Figure  4  for 

increasing  angles  of  incidence.  For  the  larger  epicentral 

distances,  we  note  a  curious  behaviour  of  the  phase  of 

the  arrival.  For  example,  for  the  epicentral  distance 

of  3  km,  which  corresponds  to  an  angle  of  incidence  of 

46.97°  for  the  ray,  the  phase  of  S.^  for  the  traces 

with  y  =-60°  and  0°  is  roughly  -tt/2  (or  3tt/2)  ,  whereas  the 
o 

phase  of  S-^  for  the  trace  with  yo=+60°  is  roughly  +tt/2. 
This  can  be  predicted  from  the  phase  curves  in  Figure  4. 

The  phase  of  S-^S.^  perfectly  elastic  case  is  roughly 

+tt/2  at  the  3  km  mark,  which  can  also  be  predicted  from  the 
phase  curve  for  the  perfectly  elastic  HlHl  coefficient. 


i  I  ■  ■  1 1  -mm  ■ 
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RAY  CODES  FOR  THE  SURFACE  POINT  SOURCE 
SYNTHETIC  SEISMOGRAMS 


Letter  Code 


Numerical  Code 


a 

b 

c 

d 

e 

f 


1  1 

1111 

111111 

11111111 

12  2  1 

1112  2  1 


129 


shown  in  Figure  14.  Hence,  if  we  assume  that  y  must  be 

o 

zero,  then  we  must  also  be  prepared  to  accept  a  phase  of 
roughly  -tt/2  for  the  anelastic  S-^S^  arrival  as  compared  to 
a  phase  of  roughly  +tt/2  for  the  elastic  arrival. 

Another  fact  pertaining  to  the  S-^S^  arrival,  obtained 
from  the  results  of  Chapter  4,  is  that  for  Yo=0°,  there 
is  no  critical  angle  for  the  transmitted  SH  ray.  As  the 
angle  of  incidence  8q  increases,  the  angle  of  transmission 
also  increases,  nearing  90°  for  0q  greater  than  about  43.29°, 
the  elastic  critical  angle,  but  never  quite  reaching  90° 
as  0q  increases  beyond  43.29°.  In  other  words,  no 
head  wave  would  be  produced  for  Yo=0°.  Critical  angles 
begin  to  exist  only  for  yo^17.43°  in  this  case. 

We  also  note  that  the  amplitude  of  the  arrival 

at  3  km  in  the  three  anelastic  cases  (Figure  13b-d)  is 
largest  for  Yq=0°,  which  again  reflects  the  behaviour  of  the 
H1H1  coefficients  in  Figure  4. 


1 . 


FIGURE  14:  HlHl  reflection  coefficient  for  the  perfectly 
elastic  case.  The  amplitude  curve  is  the  same 
as  curve  #1  in  Figure  4. 
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FOOTNOTES 


[Page  47] 

We  assume  that  jump  discontinuities  in  the  amplitude  and 
phase  curves  of  reflection/transmission  coefficients, 
such  as  in  Figure  2a,  are  unphysical  because,  refer¬ 
ring  to  the  amplitude  curve,  they  imply  a  discontinuous 
change  in  the  amplitude  of  the  reflected/transmitted 
wavefront  emanating  from  the  interface,  which  cannot  be 
accounted  for  by  any  physical  process,  since  the  assump¬ 
tion  of  homogeneity,  isotropy  and  plane  interfaces  should 
guarantee  a  continuous  wavefront.  Hence,  it  is  unaccept¬ 
able  physical  behaviour. 


[Page  89] 

Since  u(r,  t)  is  the  real  response,  we  can  replace 
( 1/2 tt )  by  (  1/tt )  Re  / o  (see  Appendix  2)  .  Also,  to 

avoid  divergence  problems  at  w  =  0 ,  we  replace  the  lower 
limit  of  the  integral  by  co i  ,  where  uj  is  small,  but 
greater  than  zero  (see  Hron  and  Kanasewich  (1971)). 
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APPENDIX  1 


FOURIER  TRANSFORM  OF  f  (t)*dg(t) 


Let  f (t)  and  g(t)  be  two  arbitrary  functions  of  time 


t.  The  Stieltjes  convolution  of  f(t)  and  g(t)  is  given  by 

t 


f(t)*dg(t)  = 


f(t-T)  d 

dx 


T 


(Al.l) 


and  its  Fourier  transform  is 


00 

f 


f  *dg  =  f(t)*dg(t)e  dt 


—  00 
00 


t 

r 


-  . ,  .  dg  (t)  ,  ,  -loot , , 

f(t-T)  — -  dx]  e  dt 

dx 


—  CO  —  co 


Changing  the  order  of  integration,  we  get 


f  *dg  = 


CO  oo 


—  co  x 


-  , .  .  -icot, ,  ,  dg  (x)  , 

f  (t-x) e  dt]  -  d 


x 


Letting  t -x=£,  the  t-integral  becomes 


x  -ioot , .  -ioox 
f ( t-x ) e  dt  =  e 


f(?)e 


X 


(Al. 2) 


(A1.3) 


(A1.4) 


Hence 


f  *dg  =  [ 


-ioot,,, 
f ( t) e  dt] 


—  ioox  dg  ( x )  , 

i  — j— -  ax 

dx 


=  [...](! 


-ICO 


-  (-ioo) 


,  ,  —  ICO X  , 

g ( x ) e  dx  } 


=  Fg 


(Al .  5 ) 


where 
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F  =  F  (03) 


g  =  g  (w) 


ICO 


-  , .  v  “icot  , . 
f  (t) e  dt 


...  -icot,^ 
g(t)  e  dt 


(A1.6) 
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APPENDIX  2 


THE  INVERSE  FOURIER  TRANSFORM 
FOR  A  REAL  RESPONSE 


Let  us  assume  that  the  response  f (t) 
receiver  is  a  real  function  of  time  t,  i.e. 


at  the 


f(t)  =  f*(t) 


where  the  *  denotes  the  complex  conjugate.  Then, 
Fourier  transform 


f  (co)  =  /  f(t)e  ^"wtdt 


we  have 


w  , 

f*  (oj)  =  /  f*  (t)elwtdt 


=  I  f  ( t)  elcotdt 


=  f  (-01) 


Hence 


f(t)  = 


1_ 

2tt 


t  -E ,  x  itotj 

J  f  (oj)  e  dw 


,  o  ^ 

=  ^  f/  ( •  •  • )  +  /  (•••)! 


OO  .  00  .  , 

=  ±_  [/  f  +  /f  (0j)ell0'da)] 

2tt  L  i 


=  L  7  [f*(w)e  lut  +  f  ( oo )  e1<iJt]  dw 

2tt 


(A2.1) 

from  its 

(A2.2) 


(A2. 3) 
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1  r  rri/  \  icot-,*  ,  rz,  .  icot-,., 

=  2^f  i  [tf(^)e  }  +  { f  (go  )  e  }  ]  dco 

o 

00 

=  ^  /  [2  Re{f  (co)  elaJt}]dw 
o 

CX) 

1  „  r  -£■ .  .  icot  , 

=  —  Re  f  (co)  e  dco 
it  1 
o 


Hence,  if  f(t)  is  real,  we  can  use  the  form  (l/7r)Re 

00 

rather  than  the  form  ( 1/ 2 tt )  /  .  The  integrand  of 

—  CO 

satisfies  (A2.3)  if  the  various  terms  are  defined 
appropriately  for  negative  frequencies. 


(A2. 4) 
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5.1) 


1  ' 


. 


